HOTFIX: Bug corrected with respect to initialization of 3D wavefields…

… + minor corrections - energy computation - screen output - Max. number of probes Signed-off-by: Guillaume Ducrozet <[email protected]>
LHEEA · Feb 23, 2015 · 0c99925 · 0c99925
1 parent 4552262
commit 0c99925
Show file tree

Hide file tree

Showing 8 changed files with 47 additions and 44 deletions.
diff --git a/sources/Benchmark/check_W_DYue.f90 b/sources/Benchmark/check_W_DYue.f90
@@ -69,7 +69,7 @@ PROGRAM check_W_DYue
 ! Adaptive Time Step Runge Kutta scheme
 TYPE(RK_parameters)        :: RK_param
 !
-REAL(RP)                   :: dt_out, dt_rk4, dt_lin, dt, dt_correc, volume, energy(4)
+REAL(RP)                   :: dt_out, dt_rk4, dt_lin, dt, dt_correc, volume, energy(3)
 REAL(RP)                   :: time_cur, time_next, h_rk, h_loc
 REAL(RP)                   :: error_old, error2
 INTEGER                    :: idx_max(2)
@@ -385,7 +385,7 @@ PROGRAM check_W_DYue
                         volume = REAL(a_eta(1,1),RP)
                         energy = calc_energy(a_eta, a_phis, da_eta)
                         !
-                        PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                        PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                         !
                         CALL fourier_2_space(a_eta, eta)
                         CALL fourier_2_space(a_phis,phis)
@@ -493,7 +493,7 @@ PROGRAM check_W_DYue
                         volume = REAL(a_eta(1,1),RP)
                         energy = calc_energy(a_eta, a_phis, da_eta)
                         !
-                        PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                        PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                         !
                         CALL fourier_2_space(a_eta, eta)
                         CALL fourier_2_space(a_phis,phis)
@@ -1211,7 +1211,7 @@ PROGRAM check_W_DYue
                                 volume = REAL(a_eta(1,1),RP)
                                 energy = calc_energy(a_eta, a_phis, da_eta)
                                 !
-                                PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                                PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                                 !
                                 CALL fourier_2_space(a_eta, eta)
                                 CALL fourier_2_space(a_phis,phis)
@@ -1318,7 +1318,7 @@ PROGRAM check_W_DYue
                                volume = REAL(a_eta(1,1),RP)
                                energy = calc_energy(a_eta, a_phis, da_eta)
                                !
-                               PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                               PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                                !
                               CALL fourier_2_space(a_eta, eta)
                               CALL fourier_2_space(a_phis,phis)
@@ -2182,7 +2182,7 @@ PROGRAM check_W_DYue
                             volume = REAL(a_eta(1,1),RP)
                             energy = calc_energy(a_eta, a_phis, da_eta)
                             !
-                            PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                            PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                             !
                             CALL fourier_2_space(a_eta, eta)
                             CALL fourier_2_space(a_phis,phis)
@@ -2291,7 +2291,7 @@ PROGRAM check_W_DYue
                             volume = REAL(a_eta(1,1),RP)
                             energy = calc_energy(a_eta, a_phis, da_eta)
                             !
-                            PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(4)
+                            PRINT*,'time cur =',time_cur,' T_stop =', T_stop_star, 'vol =', volume, 'energy = ',energy(3)
                             !
                             CALL fourier_2_space(a_eta, eta)
                             CALL fourier_2_space(a_phis,phis)
@@ -2307,13 +2307,13 @@ PROGRAM check_W_DYue
                             WRITE(77,'(A,ES9.2,A,I4,A,I4)')'ZONE T = "',time_cur,'", I=',n1o2p1,', J=',n2
                             DO i2 = n2o2p1+1, n2
                                 DO i1 = 1, n1o2p1
-                                    WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), - ky(n2-i2+2), &
+                                    WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky_n2(i2), &
                                         ABS(a_eta(i1,i2)), ABS(a_phis(i1,i2))
                                 ENDDO
                             ENDDO
                             DO i2 = 1, n2o2p1
                                 DO i1 = 1, n1o2p1
-                                    WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky(i2), &
+                                    WRITE(77,'(ES10.3,X,ES10.3,X,ES10.3,X,ES10.3)') kx(i1), ky_n2(i2), &
                                         ABS(a_eta(i1,i2)), ABS(a_phis(i1,i2))
                                 ENDDO
                             ENDDO

diff --git a/sources/HOS/HOS-ocean.f90 b/sources/HOS/HOS-ocean.f90
@@ -43,7 +43,7 @@ PROGRAM HOS_ocean
 REAL(RP)                          :: delx, dely, pioxlen, pioylen, k2, thkh
 INTEGER                           :: Mo2
 INTEGER                           :: i1, i2, j, iCPUtime, jj
-REAL(RP)                          :: energy(4)
+REAL(RP)                          :: energy(3)
 !
 REAL(RP) :: t_i, t_f, time_cur, t_tot, time_next, t_i_indiv, t_f_indiv
 REAL(RP) :: dt_out, dt_rk4, dt, h_rk, h_loc, dt_lin
@@ -284,7 +284,8 @@ PROGRAM HOS_ocean
 ! some variables useful for initialization
 ! FIXME : make it cleaner
 DO i2=1,n2o2p1
-    DO i1=1,n1o2p1
+    theta_abs(1,i2) = 0.0_rp
+    DO i1=2,n1o2p1
         k_abs(i1,i2)=SQRT(kx(i1)*kx(i1)+(ky_n2(i2)*ky_n2(i2)))
         theta_abs(i1,i2)=ATAN2(ky_n2(i2),kx(i1))
         IF (theta_abs(i1,i2) .LT. 0.0_rp) THEN
@@ -293,7 +294,8 @@ PROGRAM HOS_ocean
     ENDDO
 ENDDO
 DO i2=2,n2o2p1
-    DO i1=1,n1o2p1
+    theta_abs(1,n2-i2+1) = 0.0_rp
+    DO i1=2,n1o2p1
         k_abs(i1,n2-i2+2)=SQRT(kx(i1)*kx(i1)+(ky_n2(i2)*ky_n2(i2)))
         theta_abs(i1,n2-i2+2)=ATAN2(-ky_n2(i2),kx(i1))
         IF (theta_abs(i1,n2-i2+2) .LT. 0.0_rp) THEN
@@ -503,8 +505,8 @@ PROGRAM HOS_ocean
     !
     ! Rough estimation of peak period
     idx = MAXLOC(abs(a_eta))
-    Tp = TWOPI/(k(idx(1),idx(2))*TANH(k(idx(1),idx(2))*depth_star))
-    PRINT*,'Hs_out=',4.0_rp*SQRT(energy(4))* L_out,', T_peak=',Tp * T_out
+    Tp = TWOPI/omega_n2(idx(1),idx(2))
+    PRINT*,'Hs_out=',4.0_rp*SQRT(energy(3))* L_out,', T_peak=',Tp * T_out
     PRINT*,'***************************'
     CALL output_time_step(i_3d=i_3d, i_a=i_a_3d, i_vol=1, i_2D=i_2d, i_max=0, &
                          time=time_cur, N_stop=FLOOR(T_stop_star / dt_out), &

diff --git a/sources/HOS/energy_calc.f90 b/sources/HOS/energy_calc.f90
@@ -48,7 +48,7 @@ FUNCTION calc_energy(a_eta, a_phis, da_eta)
 ! needs     : non-dimensional gravity, g_star
 !
 COMPLEX(CP), DIMENSION(m1o2p1,m2), INTENT(IN) :: a_eta, a_phis, da_eta
-REAL(RP), DIMENSION(4)                        :: calc_energy
+REAL(RP), DIMENSION(3)                        :: calc_energy
 !
 calc_energy(:)=0.0_rp
 ! Potential energy
@@ -81,9 +81,9 @@ FUNCTION calc_energy(a_eta, a_phis, da_eta)
 !   integral
 CALL space_2_fourier_big(temp_R_Nd, temp_C_Nd)
 !
-calc_energy(3) = REAL(temp_C_Nd(1,1),RP)
-calc_energy(4) = calc_energy(1) + calc_energy(2) + calc_energy(3)
-calc_energy(1:4) = 0.5_rp * calc_energy(1:4)
+calc_energy(2) = REAL(temp_C_Nd(1,1),RP)
+calc_energy(3) = calc_energy(1) + calc_energy(2)
+calc_energy(1:3) = 0.5_rp * calc_energy(1:3)
 !
 RETURN
 !

diff --git a/sources/HOS/initial_condition.f90 b/sources/HOS/initial_condition.f90
@@ -378,15 +378,15 @@ SUBROUTINE initiate_irreg(RK_param)
 !-------------------------------------------
 !
 IMPLICIT NONE
-
+!
 REAL(RP) :: sigma, theta, E, Cj, pioxlen, pioylen
 REAL(RP) :: test, rnd, angle, angle1, angle2
-REAL(RP), DIMENSION(m1o2p1,n2) :: phi_JONSWAP
-INTEGER :: iseed,i1,i2,i_initiate_NL
-REAL(RP), DIMENSION(4)          :: energy
-TYPE(RK_parameters), INTENT(IN)         :: RK_param
+INTEGER  :: iseed,i1,i2,i_initiate_NL
+REAL(RP), DIMENSION(3)            :: energy
+TYPE(RK_parameters), INTENT(IN)   :: RK_param
 COMPLEX(CP), DIMENSION(m1o2p1,m2) :: da_eta,a_eta_temp,a_phi_temp
-REAL(RP),DIMENSION(m1o2p1,m2) :: DD_WW
+REAL(RP), DIMENSION(m1o2p1,n2) :: phi_JONSWAP
+REAL(RP),DIMENSION(m1o2p1,m2)  :: DD_WW
 REAL(RP),DIMENSION(ikp) :: ones_k
 INTEGER ::  i_dir_JSWP,jj
 REAL(RP) :: beta_min,beta_max,frr,fp_w,s_ww,Norm_DD
@@ -426,12 +426,12 @@ SUBROUTINE initiate_irreg(RK_param)
 ELSE
     DO i1=1,n1o2p1
         DO i2=1,n2
-            DD_WW(i1,i2) =1.0_rp/beta*(cos(TWOPI*theta_abs(i1,i2)/(4*beta)))**2.0_rp
+            ! With this definition of angular description, theta must be in [-pi ; pi]
+            DD_WW(i1,i2) = 1.0_rp/beta*(cos(TWOPI*ATAN2(ky_n2(i2),kx(i1))/(4*beta)))**2.0_rp
         ENDDO
     ENDDO
 ENDIF
 ! ----- End dir function
-
 iseed = 4*n1o2p1*n2_all
 !
 Cj = 3.279_rp*E_cible
@@ -567,7 +567,7 @@ SUBROUTINE initiate_irreg(RK_param)
         CALL RK_adapt_2var_3D_in_mo_lin(0, RK_param, 0.0_rp, 0.0_rp, a_phi_temp, a_eta_temp, da_eta)
         energy = calc_energy(a_eta, a_phis, da_eta)
         WRITE(*,*)'***************'
-        test = energy(4)
+        test = energy(3)
      ELSE IF(i_initiate_NL == 2)THEN
         ! FIXME : depend on wind input/dissip ?
         !CALL initiate_NL_o2
@@ -580,7 +580,7 @@ SUBROUTINE initiate_irreg(RK_param)
         ELSE
         CALL RK_adapt_2var_3D_in_mo_lin(0, RK_param, 0.0_rp, 0.0_rp, a_phi_temp, a_eta_temp, da_eta)
         energy = calc_energy(a_eta, a_phis, da_eta)
-        test = energy(4)
+        test = energy(3)
         ENDIF
         WRITE(*,*) 'E_current=', test, 'E cible =', E_cible
         E = E * E_cible /test

diff --git a/sources/HOS/output.f90 b/sources/HOS/output.f90
@@ -64,7 +64,7 @@ SUBROUTINE init_output(i_3d, i_a, i_vol, i_2D, i_max, i_prob, i_sw)
     OPEN(3,file='Results/vol_energy.dat',status='unknown')
     CALL write_input(3)
     WRITE(3,'(A)')'TITLE=" 3D volume and energy "'
-    WRITE(3,'(A)') 'VARIABLES="t", "volume", "potential", "flexural", "kinetic", "total", "dE/E_o","E_spec_dens"'
+    WRITE(3,'(A)') 'VARIABLES="t", "volume", "potential", "kinetic", "total", "dE/E_o","E_spec_dens"'
 ENDIF
 !
 IF (i_2D == 1) THEN
@@ -136,7 +136,7 @@ SUBROUTINE init_output(i_3d, i_a, i_vol, i_2D, i_max, i_prob, i_sw)
     OPEN(99,file='Results/probes.dat',status='unknown')
     CALL write_input(9)
     WRITE(99,'(A)') 'TITLE="Probes records versus time"'        !
-    WRITE(99,'(61A)') 'VARIABLES = "time" ', ('"p'//TRIM(int2str(i1))//'" ',i1=1,nprobes)
+    WRITE(99,'(101A)') 'VARIABLES = "time" ', ('"p'//TRIM(int2str(i1))//'" ',i1=1,nprobes)
 ENDIF
 !
 IF (i_sw == 1) THEN
@@ -176,7 +176,7 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
 INTEGER, INTENT(IN)               :: i_3d, i_a, i_vol, i_2D, i_max, N_stop, i_prob
 COMPLEX(CP), DIMENSION(m1o2p1,m2) :: a_eta, a_phis, da_eta
 REAL(RP), DIMENSION(m1,m2)        :: eta, phis
-REAL(RP) :: time, volume, energy(4), E_0(4), E_tot
+REAL(RP) :: time, volume, energy(3), E_0(3), E_tot
 ! Local variables
 INTEGER                        :: ii, i1, i2, it
 REAL(RP)                       :: a_1, eta_mult, min_val, max_val
@@ -286,9 +286,10 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
         ELSE
             WRITE(2,103)'ZONE T = "',time*T_out,'", I=',n1o2p1,', J=',n2
         ENDIF
+        ! Negatives ky first
         DO i2 = n2o2p1+1, n2
             DO i1 = 1, n1o2p1
-                WRITE(2,202) kx(i1)/L_out, - ky_n2(n2-i2+2)/L_out, abs_eta(i1,i2), abs_phis(i1,i2), log_eta(i1,i2), log_phis(i1,i2)
+                WRITE(2,202) kx(i1)/L_out, ky_n2(i2)/L_out, abs_eta(i1,i2), abs_phis(i1,i2), log_eta(i1,i2), log_phis(i1,i2)
             ENDDO
         ENDDO
         DO i2 = 1, n2o2p1
@@ -321,10 +322,10 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
 IF (i_vol == 1) THEN
     volume = volume * xlen_star
     IF (n2 /= 1) volume = volume * ylen_star
-    IF (ABS(E_o(4)) > tiny) THEN
-        WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,4), ABS(energy(4)-E_0(4))/E_0(4),E_tot
+    IF (ABS(E_0(3)) > tiny) THEN
+        WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,3), ABS(energy(3)-E_0(3))/E_0(3),E_tot
     ELSE
-        WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,4), 0.0_rp,E_tot
+        WRITE(3,303) time*T_out, volume*L_out, (energy(i1),i1=1,3), 0.0_rp,E_tot
     ENDIF
    303 FORMAT(7 (ES12.5,X),ES12.5)
 ENDIF
@@ -374,7 +375,11 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
         ENDDO
     ENDDO
     ! For the output of probes, use Hs_real and Tp_real...
-    WRITE(99,'(6(ES13.5,X))') time * T_out, (eta_probe(ii)* L_out, ii=1,nprobes)
+    WRITE(99,'(101(ES13.5,X))') time * T_out, (eta_probe(ii)* L_out, ii=1,nprobes)
+    IF (nprobes.GT.100) THEN
+        PRINT*, 'Change writing format in probes.dat file'
+        STOP
+    ENDIF
 ENDIF
 !
 IF (i_sw == 1) THEN ! time=0 has to be saved...

diff --git a/sources/HOS/variables_3D.f90 b/sources/HOS/variables_3D.f90
@@ -126,7 +126,7 @@ MODULE variables_3D
 REAL(RP) :: g_star, xlen_star, ylen_star, T_stop_star, f_out_star, depth_star
 !
 ! Volume and energy
-REAL(RP) :: volume, E_o(4)
+REAL(RP) :: volume, E_o(3)
 !
 ! Output numbers
 INTEGER :: i_3d, i_a_3d, i_2d, i_prob

diff --git a/sources/PostProcessing/output_post_process.f90 b/sources/PostProcessing/output_post_process.f90
@@ -1,11 +1,7 @@
 MODULE output_post_process
 !
-! This module contains the input related routines
-!  Subroutines :  read_input
-!                 read_datum
-!                 read_blank_line
-!                 build_format
-!                 error_message
+! This module contains the output related routines for initialization of output files and
+! writing at each time-step
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !

diff --git a/sources/utilities/type.f90 b/sources/utilities/type.f90
@@ -21,7 +21,7 @@ MODULE type
 !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
-! Definition of symboles for real types (RP) and complex ones (CP)
+! Definition of symbols for real types (RP) and complex ones (CP)
 !
 ! Real numbers are simple or double precision
 INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6,   37)     ! REAL32