From 188794ee96a27fd5edbeaeb178f9423991d630b9 Mon Sep 17 00:00:00 2001 From: Aron Roland Date: Sun, 6 Oct 2024 11:05:52 +0200 Subject: [PATCH] ww3_nonlinear_cg: more on the nonlinear stuff ... --- model/src/w3dispmd.F90 | 44 +++++++++++++++++++++++------------ model/src/w3profsmd_pdlib.F90 | 28 ++++++++++++++-------- model/src/w3wavemd.F90 | 5 +++- 3 files changed, 51 insertions(+), 26 deletions(-) diff --git a/model/src/w3dispmd.F90 b/model/src/w3dispmd.F90 index 6c91da560..c206e78c6 100644 --- a/model/src/w3dispmd.F90 +++ b/model/src/w3dispmd.F90 @@ -343,7 +343,7 @@ SUBROUTINE WAVNU2 (W,H,K,CG,EPS,NMAX,ICON) !/ END SUBROUTINE WAVNU2 !/ - PURE SUBROUTINE WAVNU3 (SI,H,K,CG) + SUBROUTINE WAVNU3 (SI,H,K,CG) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -452,7 +452,7 @@ PURE SUBROUTINE WAVNU3 (SI,H,K,CG) !/ END SUBROUTINE WAVNU3 - PURE SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) + SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -521,6 +521,7 @@ PURE SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) !/ ------------------------------------------------------------------- / !/ USE CONSTANTS, ONLY : GRAV, PI, TPI + USE W3GDATMD, ONLY : DMIN !!/S USE W3SERVMD, ONLY: STRACE ! IMPLICIT NONE @@ -537,8 +538,8 @@ PURE SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) INTEGER :: I1, I2 !!/S INTEGER, SAVE :: IENT = 0 REAL :: TMP, TP, CP, L, T, KD - REAL :: L0, K0, KD0, L1, TMP2, L2, TMP3, HS, HMONO - REAL :: COTH, COTH2, TANHKD, KD0NU, VBAR1, VBAR2, VBAR, UBAR, U, V + REAL :: L0, K0, KD0, L1, TMP2, L2, TMP3, HS, HMONO, K1, CG1 + REAL :: COTH, COTH2, TANHKD, KD0NU, VBAR1, VBAR2, VBAR, UBAR, U, V, DEPTH REAL, PARAMETER :: BETA1 = 1.55 REAL, PARAMETER :: BETA2 = 1.3 REAL, PARAMETER :: BETA3 = 0.216 @@ -554,27 +555,35 @@ PURE SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) !!/S CALL STRACE (IENT, 'WAVNU1') ! + DEPTH = MAX ( DMIN , DEP) + !WRITE(*,*) 'SIG, DEP', SIG,DEPTH + TP = SIG/ZPI - KD0 = ZPI*ZPI*DEP/GRAV*TP*TP + KD0 = ZPI*ZPI*DEPTH/GRAV*TP*TP TMP = 1.55 + 1.3*KD0 + 0.216*KD0*KD0 KD = KD0 * (1 + KD0**1.09 * 1./EXP(MIN(KDMAX,TMP))) / SQRT(TANH(MIN(KDMAX,KD0))) - K0 = KD/DEP - CG = 0.5*(1+(2*KD/SINH(MIN(KDMAX,2*KD))))*SIG/K0 + K1 = KD/DEPTH + CG1 = 0.5*(1+(2*KD/SINH(MIN(KDMAX,2*KD))))*SIG/K1 - HS = SQRT(4*ALOCAL*TPI*SIG/CG) + HS = SQRT(4*ALOCAL*TPI*SIG/CG1) HMONO = SQRT(2.)*HS + IF (HS .GT. 0.) WRITE(*,*) HS, HMONO, ALOCAL + T = ZPI/SIG L0 = GRAV*T*T/ZPI K0 = ZPI/L0 - KD0 = K0*DEP + KD0 = K0*DEPTH + !WRITE(*,*) 'HS, HMONO, T, L0, K0, KD0', HS, HMONO, T, L0, K0, KD0, CG - L1 = ZPI/K + L1 = ZPI/K1 TMP = KD0**(0.5*NU) TMP2 = 1.0/TANH(MIN(KDMAX,TMP)) L2 = L1 / (1 - ALPHA * (HMONO/L0) * TMP2**(2.0/NU)) K = ZPI/L2 - KD = K*DEP + !WRITE(*,'(A10,10F15.4)') 'L1, TMP, TMP2, L2, K0, K, HMONO', L1, TMP, TMP2, L2, K0, K, HMONO + + KD = K*DEPTH TMP = KD0**(0.50*NU) COTH = 1.0/TANH(TMP) COTH2 = COTH**(2.0/NU) @@ -582,29 +591,34 @@ PURE SUBROUTINE WAVNU4 (ALOCAL,SIG,DEP,K,CG) KD0NU = KD0**(0.5*NU) TMP = (-TANHKD*KD0**(0.5*NU)*COTH**2+KD0*(TANHKD-1)*(TANHKD+1)*COTH+TANHKD*(KD0)**(0.5*NU))*ALPHA*K0*HMONO*COTH**(2./NU) TMP2 = - ZPI * COTH * (-TANHKD-KD0+KD0*TANHKD**2) + !WRITE(*,'(A40,10F15.4)') 'COTH, COTH2, TANHKD, KD0NU, TMP, TMP2', COTH, COTH2, TANHKD, KD0NU, TMP, TMP2 + VBAR1 = GRAV * PI * (TMP+TMP2) - TMP = 1.d0/ZPI*ALPHA*K0*DEP*COTH2 + TMP = 1.d0/ZPI*ALPHA*K0*HMONO*COTH2 TMP2 = SQRT(GRAV*K0*TANH(KD0)/(1.d0-TMP)) TMP3 = (-ZPI+ALPHA*K0*HS*COTH2)**2*COTH VBAR2 = 1.d0/(TMP2*TMP3) + !WRITE(*,'(A40,10F15.4)') 'VBAR1, TMP, TMP2, TMP3, VBAR2', VBAR1, TMP, TMP2, TMP3, VBAR2 + U = SQRT(GRAV * K * tanh(KD) ) - UBAR = 0.5 * (GRAV * tanh(KD) + GRAV * K * (1 - tanh(KD)**2)*DEP) / U + UBAR = 0.5 * (GRAV * tanh(KD) + GRAV * K * (1 - tanh(KD)**2)*DEPTH) / U V = SQRT(1.d0/(1.d0-TMP)) TMP2 = (-0.5 * SQRT(2.d0) * PI * ALPHA * HS * COTH2) * (-COTH-KD0**(0.5*NU)+KD0**(0.5*NU)*COTH**2) TMP3 = 1.d0/(SQRT(PI/(ZPI-ALPHA*K0*HMONO*COTH2))*(-ZPI+ALPHA*K0*HMONO*COTH2)**2*COTH) VBAR = TMP2*TMP3 CG = U * VBAR + V * UBAR + !WRITE(*,'(A10,10F15.4)') 'FINAL RESULT', U, UBAR, V, VBAR, K, CG ! RETURN !/ - !/ End of WAVNU3 ----------------------------------------------------- / + !/ End of WAVNU4 ----------------------------------------------------- / !/ END SUBROUTINE WAVNU4 - PURE SUBROUTINE WAVNU_LOCAL (SIG,DW,WNL,CGL) + SUBROUTINE WAVNU_LOCAL (SIG,DW,WNL,CGL) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | diff --git a/model/src/w3profsmd_pdlib.F90 b/model/src/w3profsmd_pdlib.F90 index e70a414de..20cbd2077 100644 --- a/model/src/w3profsmd_pdlib.F90 +++ b/model/src/w3profsmd_pdlib.F90 @@ -2732,7 +2732,7 @@ SUBROUTINE CHECK_ARRAY_INTEGRAL_NX_R8(TheARR, string, maxidx) END SUBROUTINE CHECK_ARRAY_INTEGRAL_NX_R8 !/ ------------------------------------------------------------------- / !/ ------------------------------------------------------------------- / - SUBROUTINE PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC ) + SUBROUTINE PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC, ITSUB ) !/ ------------------------------------------------------------------- / !/ !/ +-----------------------------------+ @@ -2786,14 +2786,14 @@ SUBROUTINE PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC ) USE W3GDATMD, only: B_JGS_USE_JACOBI LOGICAL, INTENT(IN) :: LCALC - INTEGER, INTENT(IN) :: IMOD + INTEGER, INTENT(IN) :: IMOD, ITSUB REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY #ifdef W3_DEBUGSOLVER WRITE(740+IAPROC,*) 'B_JGS_USE_JACOBI=', B_JGS_USE_JACOBI FLUSH(740+IAPROC) #endif IF (B_JGS_USE_JACOBI) THEN - CALL PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) + CALL PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC, ITSUB) RETURN END IF WRITE(*,*) 'Error: You need to use with JGS_USE_JACOBI' @@ -3523,7 +3523,7 @@ SUBROUTINE calcARRAY_JACOBI(DTG,FACX,FACY,VGX,VGY) !/ END SUBROUTINE calcARRAY_JACOBI !/ ------------------------------------------------------------------- / - SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) + SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY,ITSUB) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -3618,6 +3618,7 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) USE W3STR1MD #endif REAL, INTENT(in) :: DTG, FACX, FACY, VGX, VGY + INTEGER, INTENT(IN) :: ITSUB INTEGER :: IP, ISP, ISEA, IP_glob INTEGER :: idx, IS INTEGER :: I, J, ITH, IK, J2 @@ -3669,10 +3670,17 @@ SUBROUTINE calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) IP_GLOB = IPLG(IP) !#ifdef NOCGTABLE - !CALL WAVNU_LOCAL(SIG(IK),DW(IP_GLOB),WN1,CG1) - CALL WAVNU4 (VA(ISP,IP),SIG(IK),DW(IP_GLOB),WN1,CG1) + IF (ITSUB .LT. 20) THEN + CALL WAVNU_LOCAL(SIG(IK),DW(IP_GLOB),WN1,CG1) + !WRITE(*,*) 'LINEAR', WN1, CG1 + ELSE + CALL WAVNU4 (VA(ISP,IP),SIG(IK),DW(IP_GLOB),WN1,CG1) + !WRITE(*,*) 'NONLINEAR', WN1, CG1 + ENDIF + + !WRITE(*,*) VA(ISP,IP),SIG(IK),DW(IP_GLOB),WN1,CG1 !#else -! CG1 = CG(IK,IP_GLOB) + !CG1 = CG(IK,IP_GLOB) !#endif CXY(1,IP) = CCOS * CG1/CLATS(IP_GLOB) CXY(2,IP) = CSIN * CG1 @@ -5455,7 +5463,7 @@ SUBROUTINE ACTION_LIMITER_LOCAL(IP,ACLOC,ACOLD, DTG) ENDIF END SUBROUTINE ACTION_LIMITER_LOCAL !/ ------------------------------------------------------------------- / - SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC) + SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC, ITSUB) !/ !/ +-----------------------------------+ !/ | WAVEWATCH III NOAA/NCEP | @@ -5551,7 +5559,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL #endif implicit none LOGICAL, INTENT(IN) :: LCALC - INTEGER, INTENT(IN) :: IMOD + INTEGER, INTENT(IN) :: IMOD, ITSUB REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY ! INTEGER :: IP, ISP, ITH, IK, JSEA, ISEA, IP_glob, IS0 @@ -5727,7 +5735,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL ! geographical advection ! IF (IMEM == 1) THEN - call calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY) + call calcARRAY_JACOBI_VEC(DTG,FACX,FACY,VGX,VGY,ITSUB) ENDIF #ifdef W3_DEBUGSOLVER diff --git a/model/src/w3wavemd.F90 b/model/src/w3wavemd.F90 index ce4ac4479..c8d9adcc7 100644 --- a/model/src/w3wavemd.F90 +++ b/model/src/w3wavemd.F90 @@ -1035,6 +1035,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & ! DO IT = IT0, NT + #ifdef W3_TIMINGS CALL PRINT_MY_TIME("Begin of IT loop") #endif @@ -1060,6 +1061,8 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 0') ! ITIME = ITIME + 1 + + WRITE(*,*) 'ITIME = ', ITIME ! DTG = REAL(NINT(DTGA+DTRES+0.0001)) DTRES = DTRES + DTGA - DTG @@ -1843,7 +1846,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT & CALL ALL_VA_INTEGRAL_PRINT(IMOD, "Before Block implicit", 1) #endif #ifdef W3_PDLIB - CALL PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE ) + CALL PDLIB_W3XYPUG_BLOCK_IMPLICIT(IMOD, FACX, FACX, DTG, VGX, VGY, UGDTUPDATE, ITIME ) IF (IAPROC == 1) WRITE(*,*) 'SUM B_JAC W3WAVE AFTER PDLIB_W3XYPUG_BLOCK_IMPLICIT', SUM(B_JAC) #endif #ifdef W3_PDLIB