Fix failure of restart reproducability in cases with variable water level and deep bathymetry. In subroutine W3ULEV KDCHK computed from bathymetric depth rather than total water level. (NOAA-EMC#1591)

Author: kestonsmith-noaa

model/src/w3initmd.F90

@@ -514,6 +514,9 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
 #ifdef W3_DIST
     CHARACTER(LEN=12)       :: FORMAT
 #endif
+    REAL                    :: DEPTHbat,WNbat,CGbat,KDCHCK
+    REAL                    :: KDMAX = 4.
+
     CHARACTER(LEN=23)       :: DTME21
     CHARACTER(LEN=30)       :: LFILE, TFILE
     integer                 :: memunit
@@ -1399,7 +1402,27 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
       WRITE (NDST,9051) IS, DEPTH
 #endif
       !
-      DO IK=0, NK+1
+      ! In shallow water, KDCHCK < KDMAX, assign time varying WN and CG as water level changes.
+      ! In deep water, KDCHCK >= KDMAX, assign time constant WN and CG based on bathymetric depth (-ZB).
+      !
+      ! Calculate the wavenumber for lowest frequency based on bathymetric depth
+      IF (IS.GT.0) THEN
+        DEPTHbat = MAX(DMIN,-ZB(IS))
+      ELSE
+        DEPTHbat = DMIN
+      END IF
+#ifdef W3_PDLIB
+      CALL WAVNU3(SIG(1),DEPTHbat,WNbat,CGbat)
+#else
+      CALL WAVNU1(SIG(1),DEPTHbat,WNbat,CGbat)
+#endif
+      KDCHCK = WNbat * DEPTHbat
+      IF ( KDCHCK .GE. KDMAX ) THEN
+        ! Use time constant CG and WN based on bathymetry for this point rather than time varying total depth.
+        DEPTH=DEPTHbat
+      ENDIF
+
+      DO IK=0, NK+1 
         !
         !         Calculate wavenumbers and group velocities.
 #ifdef W3_PDLIB

model/src/w3updtmd.F90

@@ -2201,6 +2201,7 @@ SUBROUTINE W3ULEV ( A, VA )
          OWN(NK), DWN(NK)
     REAL                    :: KDMAX = 4., RDKMIN = 0.05
     REAL                    :: WLVeff
+    REAL                    :: DEPTHbat, CGbat, WNbat
 #ifdef W3_T3
     REAL                    :: OUT(NK,NTH)
 #endif
@@ -2345,7 +2346,18 @@ SUBROUTINE W3ULEV ( A, VA )
       !
       ! 2.a Check if deep water
       !
-      KDCHCK = WN(1,ISEA) * MIN( DWO(ISEA) , DW(ISEA) )
+      !
+      ! In shallow water, KDCHCK < KDMAX, assign time varying WN and CG as water level changes.
+      ! In deep water, KDCHCK >= KDMAX, assign time constant WN and CG based on bathymetric depth (-ZB).
+      !
+      ! Calculate the wavenumber for lowest frequency based on bathymetric depth
+      DEPTHbat=MAX(DMIN,-ZB(ISEA))
+#ifdef W3_PDLIB
+      CALL WAVNU3(SIG(1),DEPTHbat,WNbat,CGbat)
+#else
+      CALL WAVNU1(SIG(1),DEPTHbat,WNbat,CGbat) 
+#endif
+      KDCHCK = WNbat * DEPTHbat
       IF ( KDCHCK .LT. KDMAX ) THEN
         !
         ! 2.b Update grid and save old grid