Work related to the down-scaling capabilities of the unstructured gri…

…d approach. (#888)


*     Consolidation and performance improvements of wave setup computation on structured and unstructured grids.
*     Performance and memory improvement of the unstructured framework.
*     Introduction of a new limiter (mixture of Hersbach & Janssen + Komen etal.)
*     Alternative computations of the group velocity for unstructured grids
*     Introduce under-relaxation for triads and wave breaking
*     Improve shallow-water source term integration for explicit source terms within the unstructured down-scaling approach.
*     Bug fixes for wetting & drying for both implicit and explicit schemes.
*     Bug fix for the wave triad interaction.
*     Introduce coastline reflection for the implicit scheme, and homogenization within the new wetting & drying scheme.
*     Introduction of Block Explicit Solver as an option for unstructured meshes.
*     Consolidation of the implicit solver in terms of memory usage, coherence.
*     debug output and further work on the Jacobi solver towards CPU coherency of the Jacobi solver.
*     Fix bugs and clean & test 2nd order time-space LAX-FCT-CRD scheme.

Co-authored-by: aliabdolali <[email protected]>
Co-authored-by: Ty Hesser <[email protected]>

model/bin/ww3_from_ftp.sh

@@ -102,9 +102,11 @@ cp -r data_regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM regtests/ww
 cp -r data_regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM regtests/ww3_tp2.14/input/toy/toy_coupled_field.nc.OASACM6
 cp -r data_regtests/ww3_tp2.14/input/toy/*.nc regtests/ww3_tp2.14/input/toy/
 cp -r data_regtests/ww3_tp2.17/input/*     regtests/ww3_tp2.17/input/
+cp -r data_regtests/ww3_tp2.19/input/*     regtests/ww3_tp2.19/input/
 cp -r data_regtests/ww3_tp2.21/input/*     regtests/ww3_tp2.21/input/
 cp -r data_regtests/mww3_test_09/input/*   regtests/mww3_test_09/input/
 cp -r data_regtests/ww3_ufs1.1/input/*     regtests/ww3_ufs1.1/input/
+cp -r data_regtests/ww3_ufs1.1/input_unstr/*     regtests/ww3_ufs1.1/input_unstr/
 cp -r data_regtests/ww3_ufs1.1/input/*.nc  regtests/ww3_ufs1.2/input/
 cp -r data_regtests/ww3_ufs1.2/input/*     regtests/ww3_ufs1.2/input/
 cp -r data_regtests/ww3_ufs1.3/input/*nc   regtests/ww3_ufs1.3/input/

model/inp/ww3_grid.inp

@@ -305,30 +305,34 @@ $   limitation and the GSE alleviation.
 $
 $ Unstructured grids ------------------------------------------------ $
 $   UNST parameters       : Namelist UNST
-$                           UGOBCAUTO : TRUE: OBC points are taken from  type 15 elements
+$                           UGOBCAUTO  : TRUE: OBC points are taken from  type 15 elements
 $                                       FALSE: OBC points must be listed in ww3_grid.inp
-$                           UGOBCDEPTH: Threshold ( < 0) depth for OBC points if UGOBCAUTO is TRUE
-$                           EXPFSN    : Activation of N scheme
+$                           UGOBCDEPTH : Threshold ( < 0) depth for OBC points if UGOBCAUTO is TRUE
+$                           UGOBCFILE : Reading boundary files from a file
+$                           EXPFSN  : Activation of N scheme (only one of the below 4, True - Active, False - not active)
 $                           EXPFSPSI  : Activation of PSI scheme
 $                           EXPFSFCT  : Activation of FCT scheme
-$                           IMPFSN    : Activation of N implicit scheme
+$                           IMPFSN  : Activation of N implicit scheme
+$                           EXPTOTAL  : Activation of the Block explicit N scheme solver
 $                           IMPTOTAL  : Activation of fully implicit scheme | Non splitting
-$                           EXPTOTAL  : Turn on implicit refraction (only with imptotal)
-$                           IMPREFRACTION  : Turn on implicit freq. shift (only with imptotal)
+$                           IMPREFRACTION : Turn on implicit freq. shift (only with imptotal)
 $                           IMPFREQSHIFT  : Turn on implicit freq. shift terms (only with imptotal)
-$                           IMPSOURCE  : Turn on implicit source terms (only with imptotal)
+$                           IMPSOURCE     : Turn on implicit source terms (only with imptotal)
 $                           JGS_TERMINATE_MAXITER  : max. Number of iterations
-$                           JGS_TERMINATE_DIFFERENCE  : terminate based on the total change of wave action
-$                           JGS_TERMINATE_NORM  : terminate based on the norm of the solution
-$                           JGS_USE_JACOBI  : Use Jacobi solver for imptotal
-$                           JGS_BLOCK_GAUSS_SEIDEL  : Use Block Gauss Seidel method for imptotal
+$                           JGS_TERMINATE_DIFFERENCE  : Terminate based on the total change of the unweightet sum of wave action
+$                           JGS_TERMINATE_NORM  : Terminate based on the norm of the solution
+$                           JGS_USE_JACOBI  : Use Jacobi solver family
+$                           JGS_BLOCK_GAUSS_SEIDEL  : Use Block Gauss Seidel method for imptotal instead of the conservative jacobi iterator.
 $                           JGS_MAXITER  : max. Number of solver iterations
 $                           JGS_PMIN  : % of grid points that do not need to converge during solver iteration.
 $                           JGS_DIFF_THR  : implicit solver threshold for JGS_TERMINATE_DIFFERENCE
 $                           JGS_NORM_THR  : terminate based on the norm of the solution
+$                           JGS_LIMITER  : use total (quasi-steady: limits whole equation) instead of local limiter (un-steady: limits only source terms)
+$                           JGS_LIMITER_FUNC  : 1 - old limiter; 2 - alternatnive limiter
 $                           SETUP_APPLY_WLV  : Compute wave setup (experimental)
 $                           SOLVERTHR_SETUP  : Solver threshold for setup computations
 $                           CRIT_DEP_SETUP  : Critical depths for setup computations
+
 $
 $ SMC grid propagation    : Namelist PSMC and default values
 $                           CFLSM  : Maximum CFL no. for propagation, 0.7

model/nml/namelists.nml

@@ -224,7 +224,7 @@ $                                            w3sic4md.ftn
 $
 $ Triad nonlinear interactions - - - - - - - - - - - - - - - - - - - -
 $   Lumped Triad Interaction (LTA) : Namelist STR1 (To be implemented)
-$                           PTRIAD1 : Proportionality coefficient (default 0.05)
+$                           PTRIAD1 : Proportionality coefficient (default 1.)
 $                           PTRIAD2 : Multiple of Tm01 up to which interaction
 $                                     is computed (2.5)
 $                           PTRIAD3 : Ursell upper limit for computing
@@ -284,27 +284,30 @@ $   limitation and the GSE alleviation.
 $
 $ Unstructured grids ------------------------------------------------ $
 $   UNST parameters       : Namelist UNST
-$                           UGOBCAUTO : TRUE: OBC points are taken from  type 15 elements
+$                           UGOBCAUTO  : TRUE: OBC points are taken from  type 15 elements
 $                                       FALSE: OBC points must be listed in ww3_grid.inp
-$                           UGOBCDEPTH: Threshold ( < 0) depth for OBC points if UGOBCAUTO is TRUE
-$                           EXPFSN    : Activation of N scheme
+$                           UGOBCDEPTH : Threshold ( < 0) depth for OBC points if UGOBCAUTO is TRUE
+$                           UGOBCFILE : Reading boundary files from a file  
+$                           EXPFSN  : Activation of N scheme (only one of the below 4, True - Active, False - not active) 
 $                           EXPFSPSI  : Activation of PSI scheme
 $                           EXPFSFCT  : Activation of FCT scheme
-$                           IMPFSN    : Activation of N implicit scheme
+$                           IMPFSN  : Activation of N implicit scheme
+$                           EXPTOTAL  : Activation of the Block explicit N scheme solver 
 $                           IMPTOTAL  : Activation of fully implicit scheme | Non splitting
-$                           EXPTOTAL  : Turn on implicit refraction (only with imptotal)
-$                           IMPREFRACTION  : Turn on implicit freq. shift (only with imptotal)
+$                           IMPREFRACTION : Turn on implicit freq. shift (only with imptotal)
 $                           IMPFREQSHIFT  : Turn on implicit freq. shift terms (only with imptotal)
-$                           IMPSOURCE  : Turn on implicit source terms (only with imptotal)
+$                           IMPSOURCE     : Turn on implicit source terms (only with imptotal)
 $                           JGS_TERMINATE_MAXITER  : max. Number of iterations
-$                           JGS_TERMINATE_DIFFERENCE  : terminate based on the total change of wave action
-$                           JGS_TERMINATE_NORM  : terminate based on the norm of the solution
-$                           JGS_USE_JACOBI  : Use Jacobi solver for imptotal
-$                           JGS_BLOCK_GAUSS_SEIDEL  : Use Block Gauss Seidel method for imptotal
+$                           JGS_TERMINATE_DIFFERENCE  : Terminate based on the total change of the unweightet sum of wave action
+$                           JGS_TERMINATE_NORM  : Terminate based on the norm of the solution
+$                           JGS_USE_JACOBI  : Use Jacobi solver family 
+$                           JGS_BLOCK_GAUSS_SEIDEL  : Use Block Gauss Seidel method for imptotal instead of the conservative jacobi iterator. 
 $                           JGS_MAXITER  : max. Number of solver iterations
 $                           JGS_PMIN  : % of grid points that do not need to converge during solver iteration.
 $                           JGS_DIFF_THR  : implicit solver threshold for JGS_TERMINATE_DIFFERENCE
 $                           JGS_NORM_THR  : terminate based on the norm of the solution
+$                           JGS_LIMITER  : use total (quasi-steady: limits whole equation) instead of local limiter (un-steady: limits only source terms)
+$                           JGS_LIMITER_FUNC  : 1 - old limiter; 2 - alternatnive limiter 
 $                           SETUP_APPLY_WLV  : Compute wave setup (experimental)
 $                           SOLVERTHR_SETUP  : Solver threshold for setup computations
 $                           CRIT_DEP_SETUP  : Critical depths for setup computations

model/src/PDLIB/yowexchangeModule.F90

@@ -107,6 +107,9 @@ module yowExchangeModule
   !> number of the second dimension for exchange
   integer, public :: n2ndDim = 1
 
+  !> number of the second dimension for exchange (nth only for wave model) 
+  integer, public :: nnthDim = 1
+
 
 contains
 

model/src/PDLIB/yowpdlibmain.F90

@@ -1309,8 +1309,8 @@ subroutine ComputeTRIA_IEN_SI_CCON
     use yowNodepool,    only: x, y, z, PDLIB_SI, PDLIB_IEN, PDLIB_TRIA, PDLIB_CCON, PDLIB_TRIA03
 
     integer I1, I2, I3, stat, IE, NI(3)
-    real(rkind) :: DXP1, DXP2, DXP3, DYP1, DYP2, DYP3, DBLTMP, TRIA03
-    logical     :: CROSSES_DATELINE
+    real  :: DXP1, DXP2, DXP3, DYP1, DYP2, DYP3, DBLTMP, TRIA03
+    logical :: CROSSES_DATELINE
 
     allocate(PDLIB_SI(npa), PDLIB_CCON(npa), PDLIB_IEN(6,ne), PDLIB_TRIA(ne), PDLIB_TRIA03(ne), stat=stat)
     if(stat/=0) call parallel_abort('SI allocation failure')
@@ -1354,7 +1354,7 @@ subroutine ComputeTRIA_IEN_SI_CCON
       PDLIB_CCON(I1) = PDLIB_CCON(I1) + 1
       PDLIB_CCON(I2) = PDLIB_CCON(I2) + 1
       PDLIB_CCON(I3) = PDLIB_CCON(I3) + 1
-      TRIA03 = PDLIB_TRIA(IE)/3.d0
+      TRIA03         = PDLIB_TRIA(IE)/3.d0
       PDLIB_SI(I1) = PDLIB_SI(I1) + TRIA03
       PDLIB_SI(I2) = PDLIB_SI(I2) + TRIA03
       PDLIB_SI(I3) = PDLIB_SI(I3) + TRIA03

model/src/w3gdatmd.F90

@@ -1034,6 +1034,7 @@ MODULE W3GDATMD
     LOGICAL :: B_JGS_USE_JACOBI
     LOGICAL :: B_JGS_BLOCK_GAUSS_SEIDEL
     INTEGER :: B_JGS_MAXITER
+    INTEGER :: B_JGS_LIMITER_FUNC
     REAL*8  :: B_JGS_PMIN
     REAL*8  :: B_JGS_DIFF_THR
     REAL*8  :: B_JGS_NORM_THR
@@ -1388,6 +1389,7 @@ MODULE W3GDATMD
   LOGICAL, POINTER :: B_JGS_USE_JACOBI
   LOGICAL, POINTER :: B_JGS_BLOCK_GAUSS_SEIDEL
   INTEGER, POINTER :: B_JGS_MAXITER
+  INTEGER, POINTER :: B_JGS_LIMITER_FUNC
   REAL(8), POINTER :: B_JGS_PMIN
   REAL(8), POINTER :: B_JGS_DIFF_THR
   REAL(8), POINTER :: B_JGS_NORM_THR
@@ -2805,6 +2807,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     B_JGS_USE_JACOBI => MPARS(IMOD)%SCHMS%B_JGS_USE_JACOBI
     B_JGS_BLOCK_GAUSS_SEIDEL => MPARS(IMOD)%SCHMS%B_JGS_BLOCK_GAUSS_SEIDEL
     B_JGS_MAXITER => MPARS(IMOD)%SCHMS%B_JGS_MAXITER
+    B_JGS_LIMITER_FUNC => MPARS(IMOD)%SCHMS%B_JGS_LIMITER_FUNC
     B_JGS_PMIN => MPARS(IMOD)%SCHMS%B_JGS_PMIN
     B_JGS_DIFF_THR => MPARS(IMOD)%SCHMS%B_JGS_DIFF_THR
     B_JGS_NORM_THR => MPARS(IMOD)%SCHMS%B_JGS_NORM_THR

model/src/w3gridmd.F90

@@ -901,6 +901,7 @@ MODULE W3GRIDMD
   LOGICAL :: JGS_TERMINATE_DIFFERENCE = .TRUE.
   LOGICAL :: JGS_TERMINATE_NORM = .TRUE.
   LOGICAL :: JGS_LIMITER = .FALSE.
+  INTEGER :: JGS_LIMITER_FUNC = 1
   LOGICAL :: JGS_BLOCK_GAUSS_SEIDEL = .TRUE.
   LOGICAL :: JGS_USE_JACOBI = .TRUE.
   LOGICAL :: JGS_SOURCE_NONLINEAR = .FALSE.
@@ -1078,6 +1079,7 @@ MODULE W3GRIDMD
        JGS_TERMINATE_DIFFERENCE,                  &
        JGS_TERMINATE_NORM,                        &
        JGS_LIMITER,                               &
+       JGS_LIMITER_FUNC,                          &
        JGS_USE_JACOBI,                            &
        JGS_BLOCK_GAUSS_SEIDEL,                    &
        JGS_MAXITER,                               &
@@ -2390,13 +2392,14 @@ SUBROUTINE W3GRID()
     IMPREFRACTION = .FALSE.
     IMPFREQSHIFT = .FALSE.
     IMPSOURCE = .FALSE.
-    SETUP_APPLY_WLV = .FALSE.
-    SOLVERTHR_SETUP=1E-14
+    SETUP_APPLY_WLV = .TRUE.
+    SOLVERTHR_SETUP=1E-6
     CRIT_DEP_SETUP=0.1
     JGS_TERMINATE_MAXITER = .TRUE.
     JGS_TERMINATE_DIFFERENCE = .TRUE.
     JGS_TERMINATE_NORM = .FALSE.
     JGS_LIMITER = .FALSE.
+    JGS_LIMITER_FUNC = 1
     JGS_BLOCK_GAUSS_SEIDEL = .TRUE.
     JGS_USE_JACOBI = .TRUE.
     JGS_MAXITER=100
@@ -2413,6 +2416,7 @@ SUBROUTINE W3GRID()
     B_JGS_TERMINATE_DIFFERENCE = JGS_TERMINATE_DIFFERENCE
     B_JGS_TERMINATE_NORM = JGS_TERMINATE_NORM
     B_JGS_LIMITER = JGS_LIMITER
+    B_JGS_LIMITER_FUNC = JGS_LIMITER_FUNC
     B_JGS_BLOCK_GAUSS_SEIDEL = JGS_BLOCK_GAUSS_SEIDEL
     B_JGS_MAXITER = JGS_MAXITER
     B_JGS_PMIN = JGS_PMIN
@@ -3266,6 +3270,7 @@ SUBROUTINE W3GRID()
            JGS_TERMINATE_DIFFERENCE,                   &
            JGS_TERMINATE_NORM,                         &
            JGS_LIMITER,                                &
+           JGS_LIMITER_FUNC,                           & 
            JGS_USE_JACOBI,                             &
            JGS_BLOCK_GAUSS_SEIDEL,                     &
            JGS_MAXITER,                                &
@@ -3600,6 +3605,10 @@ SUBROUTINE W3GRID()
         FSTOTALEXP = EXPTOTAL
         PNAME2 = 'N Explicit (Fluctuation Splitting) for one exchange explicit DC HPCF '
       END SELECT
+
+      IF (FSTOTALIMP .or. FSTOTALEXP) THEN
+        LPDLIB = .TRUE. 
+      ENDIF
       !
       IF (SUM(UNSTSCHEMES).GT.1) WRITE(NDSO,1035)
       WRITE (NDSO,2951) PNAME2
@@ -3628,7 +3637,7 @@ SUBROUTINE W3GRID()
       END IF
       IF (SETUP_APPLY_WLV) THEN
         DO_CHANGE_WLV = SETUP_APPLY_WLV
-        PNAME2 = ' we change WLV'
+        PNAME2 = 'Wave setup is added to the WLV'
         WRITE (NDSO,2952) PNAME2
       END IF
       SOLVERTHR_STP = SOLVERTHR_SETUP
@@ -6571,6 +6580,7 @@ SUBROUTINE W3GRID()
          ',  JGS_TERMINATE_DIFFERENCE=', L3,                    &
          ',  JGS_TERMINATE_NORM=', L3,                          &
          ',  JGS_LIMITER=', L3,                                 &
+         ',  JGS_LIMITER_FUNC=', I3,                            &
          ',  JGS_USE_JACOBI=', L3,                              &
          ',  JGS_BLOCK_GAUSS_SEIDEL=', L3,                      &
          ',  JGS_MAXITER=', I5,                                 &

model/src/w3initmd.F90

@@ -388,7 +388,7 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
          ICLOSE, ZB, TRNX, TRNY, DMIN, DTCFL, DTMAX, &
          FLCK, NK, NTH, NSPEC, SIG, GNAME
 #ifdef W3_PDLIB
-    USE W3GDATMD, ONLY : FLCTH
+    USE W3GDATMD, ONLY : FLCTH, B_JGS_BLOCK_GAUSS_SEIDEL, B_JGS_USE_JACOBI
 #endif
     USE W3WDATMD, ONLY: TIME, TLEV, TICE, TRHO, WLV, UST, USTDIR, VA
     USE W3ODATMD, ONLY: NDSO, NDSE, NDST, SCREEN, NDS, NTPROC,      &
@@ -421,7 +421,7 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
     USE W3GDATMD, ONLY: GTYPE, UNGTYPE
 #ifdef W3_PDLIB
     USE PDLIB_W3PROFSMD, ONLY : PDLIB_MAPSTA_INIT, SET_IOBDP_PDLIB, PDLIB_IOBP_INIT, SET_IOBPA_PDLIB
-    USE PDLIB_W3PROFSMD, ONLY : BLOCK_SOLVER_INIT, PDLIB_INIT, DEALLOCATE_PDLIB_GLOBAL
+    USE PDLIB_W3PROFSMD, ONLY : BLOCK_SOLVER_INIT, BLOCK_SOLVER_EXPLICIT_INIT, PDLIB_INIT, DEALLOCATE_PDLIB_GLOBAL
     use yowDatapool, only: istatus
 #endif
 #ifdef W3_SETUP
@@ -639,12 +639,23 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
     IF (FSTOTALIMP .and. .NOT. LPDLIB) THEN
       WRITE(NDSE,*) 'IMPTOTAL is selected'
       WRITE(NDSE,*) 'But PDLIB is not'
-      STOP 'Stop, case 1'
+      CALL FLUSH(NDSE) 
+      STOP 
     ELSE IF (FSTOTALEXP .and. .NOT. LPDLIB) THEN
       WRITE(NDSE,*) 'EXPTOTAL is selected'
       WRITE(NDSE,*) 'But PDLIB is not'
-      STOP 'Stop, case 1'
+      CALL FLUSH(NDSE) 
+      STOP 
     END IF
+#ifdef W3_PDLIB
+    IF (B_JGS_BLOCK_GAUSS_SEIDEL .AND. .NOT. B_JGS_USE_JACOBI) THEN
+      WRITE(NDSE,*) 'B_JGS_BLOCK_GAUSS_SEIDEL is used but the Jacobi solver is not choosen'
+      WRITE(NDSE,*) 'Please set JGS_USE_JACOBI .eqv. .true.'
+      CALL FLUSH(NDSE) 
+      STOP 
+    ENDIF
+#endif
+
     !
     ! 1.c Open files without unpacking MDS ,,,
     !
@@ -792,17 +803,12 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
 
 #ifdef W3_PDLIB
     IF ((IAPROC .LE. NAPROC).and.(GTYPE .eq. UNGTYPE)) THEN
-#endif
-
-#ifdef W3_PDLIB
       CALL BLOCK_SOLVER_INIT(IMOD)
       CALL PDLIB_IOBP_INIT(IMOD)
       CALL SET_IOBPA_PDLIB
-#endif
-
-#ifdef W3_PDLIB
-    ELSE IF (FSTOTALEXP) THEN
-      !AR: To do here the blocksolver ...
+      IF (FSTOTALEXP) THEN
+        CALL BLOCK_SOLVER_EXPLICIT_INIT()
+      ENDIF
     ENDIF
 #endif
 
@@ -1264,6 +1270,10 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
     !Li   DO IY=1, NY
     !Li     DO IX=1, NX
     !Li       ISEA   = MAPFS(IY,IX)
+#ifdef W3_DEBUGSTP
+    max_val = 0
+    min_val = 0
+#endif
     DO ISEA=1, NSEA
       IX = MAPSF(ISEA,1)
       IY = MAPSF(ISEA,2)
@@ -1276,6 +1286,10 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
       IF (DO_CHANGE_WLV) THEN
         WLVeff=WLVeff + ZETA_SETUP(ISEA)
       END IF
+#endif
+#ifdef W3_DEBUGSTP
+      max_val = MAX(max_val, WLVeff)
+      min_val = MIN(min_val, WLVeff)
 #endif
       DW(ISEA) = MAX ( 0. , WLVeff-ZB(ISEA) )
 #ifdef W3_T
@@ -1288,19 +1302,33 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
       !Li     END IF
     END DO
     !Li   END DO
+#ifdef W3_DEBUGSTP
+      WRITE(740+IAPROC,*) 'w3initmd 1: max/min(WLVeff)=', max_val, min_val
+      FLUSH(740+IAPROC)
+      max_val = 0
+      min_val = 0
+#endif
     DO JSEA=1, NSEAL
       CALL INIT_GET_ISEA(ISEA, JSEA)
       WLVeff=WLV(ISEA)
 #ifdef W3_SETUP
       IF (DO_CHANGE_WLV) THEN
         WLVeff=WLVeff + ZETA_SETUP(ISEA)
       END IF
+#endif
+#ifdef W3_DEBUGSTP
+      max_val = MAX(max_val, WLVeff)
+      min_val = MIN(min_val, WLVeff)
 #endif
       DW(ISEA) = MAX ( 0. , WLVeff-ZB(ISEA) )
       IF ( WLVeff-ZB(ISEA) .LE.0. ) THEN
         VA(:,JSEA) = 0.
       END IF
     END DO
+#ifdef W3_DEBUGSTP
+    WRITE(740+IAPROC,*) 'w3initmd 2: max/min(WLVeff)=', max_val, min_val
+    FLUSH(740+IAPROC)
+#endif
     !
 #ifdef W3_PDLIB
     IF ( IAPROC .LE. NAPROC ) THEN
@@ -1353,8 +1381,12 @@ SUBROUTINE W3INIT ( IMOD, IsMulti, FEXT, MDS, MTRACE, ODAT, FLGRD,  FLGR2, FLGD,
       !
       DO IK=0, NK+1
         !
-        !  Calculate wavenumbers and group velocities.
+        !         Calculate wavenumbers and group velocities.
+#ifdef W3_PDLIB
+        CALL WAVNU3(SIG(IK),DEPTH,WN(IK,IS),CG(IK,IS))
+#else
         CALL WAVNU1(SIG(IK),DEPTH,WN(IK,IS),CG(IK,IS))
+#endif
         !
 #ifdef W3_T1
         WRITE (NDST,9052) IK, TPI/SIG(IK), WN(IK,IS), CG(IK,IS)

model/src/w3iogomd.F90

@@ -1533,7 +1533,7 @@ SUBROUTINE W3OUTG ( A, FLPART, FLOUTG, FLOUTG2 )
         CALL INIT_GET_ISEA(ISEA, JSEA)
         FACTOR       = DDEN(IK) / CG(IK,ISEA)
         EBD(IK,JSEA) = AB(JSEA) * FACTOR
-        ET (JSEA)    = ET (JSEA) + EBD(IK,JSEA)
+        ET(JSEA)     = ET(JSEA) + EBD(IK,JSEA)
 #ifdef W3_IG1
         IF (IK.EQ.NINT(IGPARS(5))) HSIG(JSEA) = 4*SQRT(ET(JSEA))
 #endif

model/src/w3iorsmd.F90

@@ -769,7 +769,8 @@ SUBROUTINE W3IORS ( INXOUT, NDSR, DUMFPI, IMOD, FLRSTRT )
               !              Include remainder values (switch to record format) ---- *
               JSEA = NSEAL_MIN + 1
               IF ( JSEA.EQ.NSEAL ) THEN
-                ISEA = IAPROC + (JSEA - 1) * NAPROC
+                !ISEA = IAPROC + (JSEA - 1) * NAPROC
+                CALL INIT_GET_ISEA(ISEA, JSEA)
                 NREC = ISEA + 2
                 RPOS = 1_8 + LRECL*(NREC-1_8)
                 READ (NDSR, POS=RPOS, ERR=802, IOSTAT=IERR)          &