HOTFIX: Velocity comp. for irregular waves in finite depth+minor corr…

…ections - g_star is not equal to one in irregular seas with finite depth - outputs adjusted - test to ensure that regular waves are correctly used - useless display_velocity subroutine removed Signed-off-by: Guillaume Ducrozet <guillaume.ducrozet@ec-nantes.fr>
LHEEA · Nov 30, 2015 · 6afa67b · 6afa67b
1 parent b833648
commit 6afa67b
Show file tree

Hide file tree

Showing 4 changed files with 21 additions and 30 deletions.
diff --git a/sources/HOS/HOS-ocean.f90 b/sources/HOS/HOS-ocean.f90
@@ -362,8 +362,6 @@ PROGRAM HOS_ocean
     ENDDO
 ENDDO
 !
-CALL display_velocities()
-!
 ! Initial conditions
 !
 ! Evaluates the modal amplitudes (FT)
@@ -600,35 +598,21 @@ SUBROUTINE check_slope
 !
 slopex_max = MAXVAL(ABS(etax))
 slopey_max = MAXVAL(ABS(etay))
-! WRITE(*,*) slopex_max
 !
 IF (slopex_max > threshold) THEN
-    WRITE(*,'(A,2(ES11.3,X))') 'Ca va trancher en x...', time_cur, slopex_max
+    WRITE(*,'(A,2(ES11.3,X))') 'Slope along x is too large...', time_cur, slopex_max
     STOP
 ENDIF
 !
 IF (slopey_max > threshold) THEN
-    WRITE(*,'(A)') 'Ca va trancher en y...'
+    WRITE(*,'(A)') 'Slope along y is too large...'
     STOP
 ENDIF
 !
 END SUBROUTINE check_slope
 !
 !
 !
-SUBROUTINE display_velocities()
-!
-IMPLICIT NONE
-!
-WRITE(*,'(A)') 'Parameters for the x-direction'
-IF (depth < 1.0E14_rp) THEN
-    WRITE(*,'(3X, A, ES12.5)')  'Shallow water velocity             c_0      =',SQRT(g_star*depth_star)*L/T
-ENDIF
-!
-END SUBROUTINE display_velocities
-!
-!
-!
 SUBROUTINE display_parameters()
 !
 IMPLICIT NONE

diff --git a/sources/HOS/initial_condition.f90 b/sources/HOS/initial_condition.f90
@@ -76,7 +76,7 @@ SUBROUTINE initiate_parameters()
     CASE (1,2,21)
         WRITE(*,'(A,I3)') 'initiate_parameters: No parameters to set with i_case=',i_case
     CASE(81,82,83,84,809)
-        WRITE(*,'(A)') 'initiate_parameters: Fructus case'
+        WRITE(*,'(A)') 'initiate_parameters: Rienecker & Fenton case'
         SELECT CASE (i_case)
                 CASE (81)
                 filename = 'waverf_L628_inf_ka01_N15_30.cof'
@@ -90,7 +90,11 @@ SUBROUTINE initiate_parameters()
                 filename = 'waverf_L628_inf_ka009_N50_100.cof'
         END SELECT
         ! Depth is infinite in those cases
-        depth = 1.0e15_rp
+        ! Check that correct depth is used in input file
+        IF (ABS(depth-1.0e15_rp) >= epsilon(1.0e15_rp)) THEN
+            WRITE(*,'(A)') 'For this test case, infinite water depth has to be used'
+            STOP
+        ENDIF
         !
         CALL read_RF_data(filename, RF_obj, .TRUE.)
         !

diff --git a/sources/HOS/output.f90 b/sources/HOS/output.f90
@@ -185,7 +185,7 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
 !
 eta_mult = eta_out
 IF (i_3D == 1) THEN
-    ! Analytical integration of the linear part
+    ! Transform modal amplitudes to physical variables
     CALL fourier_2_space(a_eta,  eta)
     CALL fourier_2_space(a_phis, phis)
     !
@@ -225,11 +225,14 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
         ENDDO
     ENDIF
    102 FORMAT(3(ES12.5,X),ES12.5)
-   103 FORMAT(A,F9.2,A,I5,A,I5)
+   103 FORMAT(A,ES12.5,A,I5,A,I5)
    104 FORMAT((ES12.5,X),ES12.5)
 ENDIF
 !
 IF (i_2D == 1) THEN
+    ! Transform modal amplitude to physical variable
+    CALL fourier_2_space(a_eta,  eta)
+    !
     IF (ABS(time) <= tiny) THEN
         IF (tecplot == 11) THEN
             WRITE(66,603)'ZONE SOLUTIONTIME = ',time*T_out,', I=',n1
@@ -250,10 +253,10 @@ SUBROUTINE output_time_step(i_3d, i_a, i_vol, i_2D, i_max, time, N_stop, a_eta,
         ENDDO
     ENDIF
    602 FORMAT(ES12.5,X,ES12.5)
-   603 FORMAT(A,F9.2,A,I5)
+   603 FORMAT(A,ES12.5,A,I5)
    604 FORMAT(ES12.5)
 ELSE IF (i_2D == 2) THEN
-    ! Analytical integration of the linear part
+    ! Transform modal amplitude to physical variable
     CALL fourier_2_space(a_eta,  eta)
     !
     ! Output of the free surface elevation versus space (size control of the output file)

diff --git a/sources/HOS/velocities.f90 b/sources/HOS/velocities.f90
@@ -98,7 +98,7 @@ SUBROUTINE HOSvel2(meta2,M_HOSvel,a_eta_l,a_phis_l,time_current)
 DO i2=1,n2
     DO i1=1,n1
         deta3(i1,i2)  = (deta3(i1,i2) + W1(i1,i2))
-        dphis3(i1,i2) = (dphis3(i1,i2) - eta3(i1,i2))
+        dphis3(i1,i2) = (dphis3(i1,i2) - g_star*eta3(i1,i2))
     ENDDO
 ENDDO
 !
@@ -251,10 +251,10 @@ SUBROUTINE HOSvel2(meta2,M_HOSvel,a_eta_l,a_phis_l,time_current)
         vity2ref_SL(i1,i2) = phimy2(i1,i2)
         vitz2ref_SL(i1,i2) = phimz2(i1,i2)
 
-        Press_SL(i1,i2) = - eta3(i1,i2) - phit(i1,i2) +(- 0.5_rp*(phimx2(i1,i2) + phimy2(i1,i2) + phimz2(i1,i2)))    ! DFSBC OK if phimx2.NE.phimx*phimx
+        Press_SL(i1,i2) = - g_star*eta3(i1,i2) - phit(i1,i2) +(- 0.5_rp*(phimx2(i1,i2) + phimy2(i1,i2) + phimz2(i1,i2)))    ! DFSBC OK if phimx2.NE.phimx*phimx
         CCSL(i1,i2) = deta3(i1,i2) - phimxetax(i1,i2) - phimz(i1,i2)       ! CCSL OK
         ! Approximation (no dealiasing)
-        Press_bis(i1,i2) = - eta3(i1,i2) - phit(i1,i2) - 0.5_rp*(phimx(i1,i2)**2 + phimy(i1,i2)**2 + phimz(i1,i2)**2)
+        Press_bis(i1,i2) = - g_star*eta3(i1,i2) - phit(i1,i2) - 0.5_rp*(phimx(i1,i2)**2 + phimy(i1,i2)**2 + phimz(i1,i2)**2)
         CCSL_bis(i1,i2)  = deta3(i1,i2) + (phimx(i1,i2)*etax_tmp(i1,i2)+phimy(i1,i2)*etay_tmp(i1,i2)) - phimz(i1,i2) ! KFSBC OK
     ENDDO
 ENDDO
@@ -902,7 +902,7 @@ SUBROUTINE HOSvel_direct(a_eta_l,a_phis_l,time_current)
 DO i2=1,n2
     DO i1=1,n1
         deta3(i1,i2)  = (deta3(i1,i2) + W1(i1,i2))
-        dphis3(i1,i2) = (dphis3(i1,i2) - eta3(i1,i2))
+        dphis3(i1,i2) = (dphis3(i1,i2) - g_star*eta3(i1,i2))
     ENDDO
 ENDDO
 !
@@ -987,10 +987,10 @@ SUBROUTINE HOSvel_direct(a_eta_l,a_phis_l,time_current)
         vity2ref_SL(i1,i2) = phimy2(i1,i2)
         vitz2ref_SL(i1,i2) = phimz2(i1,i2)
 
-        Press_SL(i1,i2) = - eta3(i1,i2) - phit(i1,i2) +(- 0.5_rp*(phimx2(i1,i2) + phimy2(i1,i2) + phimz2(i1,i2))) ! CDSL OK si phimx2.NE.phimx*phimx
+        Press_SL(i1,i2) = - g_star*eta3(i1,i2) - phit(i1,i2) +(- 0.5_rp*(phimx2(i1,i2) + phimy2(i1,i2) + phimz2(i1,i2))) ! CDSL OK si phimx2.NE.phimx*phimx
         CCSL(i1,i2) = deta3(i1,i2) - phimxetax(i1,i2) - phimz(i1,i2)       ! CCSL OK
         ! Approximation (no dealiasing)
-        Press_bis(i1,i2) = - eta3(i1,i2) - phit(i1,i2) - 0.5_rp*(phimx(i1,i2)**2 + phimy(i1,i2)**2 + phimz(i1,i2)**2)
+        Press_bis(i1,i2) = - g_star*eta3(i1,i2) - phit(i1,i2) - 0.5_rp*(phimx(i1,i2)**2 + phimy(i1,i2)**2 + phimz(i1,i2)**2)
         CCSL_bis(i1,i2)  = deta3(i1,i2) + (phimx(i1,i2)*etax_tmp(i1,i2)+phimy(i1,i2)*etay_tmp(i1,i2)) - phimz(i1,i2)       ! CCSL OK
     ENDDO
 ENDDO