Skip to content

Commit

Permalink
ww3_nonlinear_cg: more on the nonlinear stuff ...
Browse files Browse the repository at this point in the history
  • Loading branch information
@aronroland
aronroland committed Oct 6, 2024
1 parent 31c9017 commit 188794e
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 26 deletions.
44 changes: 29 additions & 15 deletions model/src/w3dispmd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 |
Expand Down Expand Up @@ -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 |
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -554,57 +555,70 @@ 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)
TANHKD = tanh(KD0)
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 |
Expand Down
28 changes: 18 additions & 10 deletions model/src/w3profsmd_pdlib.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
!/ ------------------------------------------------------------------- /
!/
!/ +-----------------------------------+
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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 |
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 |
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion model/src/w3wavemd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 188794e

Please sign in to comment.