Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Init bugs2 #1142

Merged
merged 15 commits into from
Dec 12, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 18 additions & 6 deletions model/src/w3profsmd_pdlib.F90
Original file line number Diff line number Diff line change
Expand Up @@ -2854,12 +2854,24 @@ SUBROUTINE PDLIB_W3XYPUG_BLOCK_EXPLICIT(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
!
USE W3ODATMD, only: IAPROC
USE W3GDATMD, only: B_JGS_USE_JACOBI
USE W3TIMEMD, only: DSEC21
USE W3ODATMD, only: TBPI0, TBPIN, FLBPI
USE W3WDATMD, only: TIME

LOGICAL, INTENT(IN) :: LCALC
INTEGER, INTENT(IN) :: IMOD
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
REAL :: RD1, RD2

IF ( FLBPI ) THEN
RD1 = DSEC21 ( TBPI0, TIME )
RD2 = DSEC21 ( TBPI0, TBPIN )
ELSE
RD1=1.
RD2=0.
END IF

CALL PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
CALL PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, RD1, RD2, DTG, VGX, VGY, LCALC)
!/
!/ End of W3XYPFSN ----------------------------------------------------- /
!/
Expand Down Expand Up @@ -6328,7 +6340,7 @@ SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCAL
#endif
END SUBROUTINE PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK
!/ ------------------------------------------------------------------- /
SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, RD10, RD20, DTG, VGX, VGY, LCALC)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
Expand Down Expand Up @@ -6402,7 +6414,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)

INTEGER, INTENT(IN) :: IMOD

REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
REAL, INTENT(IN) :: FACX, FACY, DTG, VGX, VGY, RD10, RD20

REAL :: KTMP(3), UTILDE(NTH), ST(NTH,NPA)
REAL :: FL11(NTH), FL12(NTH), FL21(NTH), FL22(NTH), FL31(NTH), FL32(NTH), KKSUM(NTH,NPA)
Expand All @@ -6411,7 +6423,7 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
REAL :: KSIG(NPA), CGSIG(NPA), CXX(NTH,NPA), CYY(NTH,NPA)
REAL :: LAMBDAX(NTH), LAMBDAY(NTH)
REAL :: DTMAX(NTH), DTMAXEXP(NTH), DTMAXOUT, DTMAXGL
REAL :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2, RD10, RD20
REAL :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2
REAL :: UOLD(NTH,NPA), U(NTH,NPA)

REAL, PARAMETER :: ONESIXTH = 1.0/6.0
Expand Down Expand Up @@ -6570,8 +6582,8 @@ SUBROUTINE PDLIB_EXPLICIT_BLOCK(IMOD, FACX, FACY, DTG, VGX, VGY, LCALC)
IF ( FLBPI ) THEN
DO ITH = 1, NTH
ISP = ITH + (IK-1) * NTH
RD1 = RD10 - DTG * REAL(ITER(IK)-IT)/REAL(ITER(IK))
RD2 = RD20
RD1=RD10 - DTMAXGL * REAL(ITER(IK)-IT)/REAL(ITER(IK))
RD2=RD20
IF ( RD2 .GT. 0.001 ) THEN
RD2 = MIN(1.,MAX(0.,RD1/RD2))
RD1 = 1. - RD2
Expand Down
12 changes: 5 additions & 7 deletions model/src/w3sdb1md.F90
Original file line number Diff line number Diff line change
Expand Up @@ -187,6 +187,7 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
USE W3ODATMD, ONLY: NDST
USE W3GDATMD, ONLY: SIG
USE W3ODATMD, only : IAPROC
USE W3PARALL, only : THR
#ifdef W3_S
USE W3SERVMD, ONLY: STRACE
#endif
Expand Down Expand Up @@ -218,7 +219,7 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
INTEGER, SAVE :: IENT = 0
#endif
REAL*8 :: HM, BB, ARG, Q0, QB, B, CBJ, HRMS, EB(NK)
REAL*8 :: AUX, CBJ2, RATIO, S0, S1, THR, BR1, BR2, FAK
REAL*8 :: AUX, CBJ2, RATIO, S0, S1, BR1, BR2, FAK
REAL :: ETOT, FMEAN2
#ifdef W3_T0
REAL :: DOUT(NK,NTH)
Expand All @@ -231,12 +232,9 @@ SUBROUTINE W3SDB1 (IX, A, DEPTH, EMEAN, FMEAN, WNMEAN, CG, LBREAK, S, D )
#endif
!
! 0. Initialzations ------------------------------------------------- /
! Never touch this 4 lines below ... otherwise my exceptionhandling will not work.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was doing some testing below - and saw this "Never touch these four lines" comment. I'm assuming @aronroland that the exception handling still works.... but just thought I'd draw attention to this comment.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @JessicaMeixner-NOAA sorry for my late reply, yes this was a leftover from long ago so I cleaned the code thanks for looking at this.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @JessicaMeixner-NOAA, back to your comment on the UFS and the DB1 issue. I have check more the code and the reason is in this case not DB1 but rather my changes in w3wave, where we changed how the source terms are called in the implicit schemes for the sake of neumann boundary but also wetting and drying. I will expand a bit the text so it becomes more clear to everybody.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for replying @aronroland - and thanks for cleaning this up and all the bug fixes!

S = 0.
D = 0.

THR = DBLE(1.E-15)
IF (SUM(A) .LT. THR) RETURN
IF (EMEAN .LT. TINY(1.d0)) THEN
RETURN
ENDIF

IWB = 1
!
Expand Down
17 changes: 11 additions & 6 deletions model/src/w3srcemd.F90
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
!> @file

!> @brief Source term integration routine.
!>
!> @author H. L. Tolman
Expand Down Expand Up @@ -1244,7 +1244,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
IF (.NOT. FSSOURCE .or. LSLOC) THEN
#endif
#ifdef W3_TR1
CALL W3STR1 ( SPEC, SPECOLD, CG1, WN1, DEPTH, IX, VSTR, VDTR )
CALL W3STR1 ( SPEC, CG1, WN1, DEPTH, IX, VSTR, VDTR )
#endif
#ifdef W3_PDLIB
ENDIF
Expand Down Expand Up @@ -1534,8 +1534,13 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
DVS = SIGN(MIN(MAXDAC,ABS(DVS)),DVS)
ENDIF
PreVS = DVS / FAKS
eVS = PreVS / CG1(IK) * CLATSL
eVD = MIN(0.,VD(ISP))
IF (IOBP_LOC(JSEA) .EQ. 3) THEN
eVS = 0
eVD = 0
ELSE
eVS = PreVS / CG1(IK) * CLATSL
eVD = MIN(0.,VD(ISP))
ENDIF
B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * (eVS - eVD*SPEC(ISP)*JAC)
ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
#ifdef W3_DB1
Expand All @@ -1548,9 +1553,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
evS = -evS
evD = 2*evD
ENDIF
#endif
B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * eVS
ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
#endif

#ifdef W3_TR1
eVS = VSTR(ISP) * JAC
Expand All @@ -1562,9 +1567,9 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, &
evS = -evS
evD = 2*evD
ENDIF
#endif
B_JAC(ISP,JSEA) = B_JAC(ISP,JSEA) + SIDT * eVS
ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) = ASPAR_JAC(ISP,PDLIB_I_DIAG(JSEA)) - SIDT * eVD
#endif
END DO
END DO

Expand Down
11 changes: 6 additions & 5 deletions model/src/w3str1md.F90
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ MODULE W3STR1MD
!>
!> @author A. J. van der Westhuysen @date 13-Jan-2013
!>
SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
SUBROUTINE W3STR1 (A, CG, WN, DEPTH, IX, S, D)
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III NOAA/NCEP |
Expand Down Expand Up @@ -259,7 +259,6 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
! CG R.A. I Group velocities.
! WN R.A. I Wavenumbers.
! DEPTH Real I Mean water depth.
! EMEAN Real I Mean wave energy.
! FMEAN Real I Mean wave frequency.
! S R.A. O Source term (1-D version).
! D R.A. O Diagonal term of derivative (1-D version).
Expand Down Expand Up @@ -320,7 +319,7 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
REAL, INTENT(IN) :: CG(NK), WN(NK), DEPTH, A(NSPEC), AOLD(NSPEC)
REAL, INTENT(IN) :: CG(NK), WN(NK), DEPTH, A(NSPEC)
INTEGER, INTENT(IN) :: IX
REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC)
!/
Expand Down Expand Up @@ -391,11 +390,13 @@ SUBROUTINE W3STR1 (A, AOLD, CG, WN, DEPTH, IX, S, D)
#ifdef W3_S
CALL STRACE (IENT, 'W3STR1')
#endif

!AR: todo: check all PRX routines for differences, check original thesis of elderberky.
!
! 1. Integral over directions
!
IF (MAXVAL(A) .LT. TINY(1.)) THEN
RETURN
ENDIF

SIGM01 = 0.
EMEAN = 0.
JACEPS = 1E-12
Expand Down
10 changes: 10 additions & 0 deletions model/src/w3wavemd.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1453,6 +1453,12 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &
call print_memcheck(memunit, 'memcheck_____:'//' WW3_WAVE TIME LOOP 13')
!
#ifdef W3_PDLIB

IF (LPDLIB .and. .not. FLSOU .and. .not. FSSOURCE) THEN
B_JAC = 0.
ASPAR_JAC = 0.
ENDIF

IF (LPDLIB .and. FLSOU .and. FSSOURCE) THEN
#endif

Expand Down Expand Up @@ -1484,6 +1490,8 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &

CALL INIT_GET_ISEA(ISEA, JSEA)

IF ((IOBP_LOC(JSEA).eq.1..or.IOBP_LOC(JSEA).eq. 3).and.IOBDP_LOC(JSEA).eq.1.and.IOBPA_LOC(JSEA).eq.0) THEN

IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
DELA=1.
Expand Down Expand Up @@ -1556,6 +1564,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &
WRITE(740+IAPROC,*) ' SHAVETOT=', SHAVETOT(JSEA)
FLUSH(740+IAPROC)
#endif
ENDIF
END DO ! JSEA
END IF ! PDLIB
#endif
Expand Down Expand Up @@ -2158,6 +2167,7 @@ SUBROUTINE W3WAVE ( IMOD, ODAT, TEND, STAMP, NO_OUT &
!
DO JSEA=1, NSEAL
CALL INIT_GET_ISEA(ISEA, JSEA)

IX = MAPSF(ISEA,1)
IY = MAPSF(ISEA,2)
DELA=1.
Expand Down
Loading