Merge branch 'NOAA-EMC:develop' into bug/points180360

model/src/w3profsmd_pdlib.F90

@@ -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 ----------------------------------------------------- /
     !/
@@ -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 |
@@ -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)
@@ -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
@@ -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

model/src/w3sdb1md.F90

@@ -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
@@ -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)
@@ -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.
-    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
     !

model/src/w3sic4md.F90

@@ -89,6 +89,7 @@ MODULE W3SIC4MD
   !        *** Rogers et al. tech. rep. 2021 (RYW2021)
   !        *** Yu et al. CRST 2022
   !        *** Yu JMSE 2022
+  !        *** Meylan et al. Ocean Modeling 2021
   !
   !  6. Switches :
   !
@@ -138,6 +139,7 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
     !/    11-Jan-2024 : Method 8 added (Meylan et al. 2018)   (E. Rogers)
     !/    11-Jan-2024 : Method 9 added (Rogers et al., 2021)
     !/                                      denoted "RYW2021" (E. Rogers)
+    !/    14-Aug-2024 : Method 10 added (Meylan et al. 2021)  (E. Thomas)
     !/
     !/        FIXME   : Move field input to W3SRCE and provide
     !/     (S.Zieger)   input parameter to W3SIC1 to make the subroutine
@@ -307,6 +309,8 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
     !        suggested default is marked with "(*SD*)", for consistency
     !          with SWAN (v41.31AB or later)
     !
+    !     10) Meylan et al. 2021 (Ocean Modeling): ocean-wave attenuation 
+    !         due to scattering by sea ice floes. 
     !     ------------------------------------------------------------------
     !
     !     For all methods, the user can specify namelist
@@ -450,6 +454,8 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
     REAL, ALLOCATABLE       :: FREQ(:) ! wave frequency
     REAL, ALLOCATABLE       :: MARG1(:), MARG2(:) ! Arguments for M2
     REAL, ALLOCATABLE       :: KARG1(:), KARG2(:), KARG3(:) !Arguments for M3
+    REAL                    :: x1,x2,x3,x1sqr,x2sqr,x3sqr   !Arguments for M10
+    REAL                    :: perfour,amhb,bmhb            !Arguments for M10
     LOGICAL                 :: NML_INPUT ! if using namelist input for M2
 
     !/
@@ -699,6 +705,43 @@ SUBROUTINE W3SIC4 (A, DEPTH, CG, IX, IY, S, D)
       DO IK=1,NK
         WN_I(IK)  = Chf*(hice**mpow)*(FREQ(IK)**npow)
       END DO
+
+    CASE (10)
+      ! Cubic fit to Meylan, Horvat & Bitz 2021
+      ! ICECOEF1 is thickness
+      ! ICECOEF5 is floe size
+      ! TPI/SIG is period
+      x3=min(ICECOEF1,3.5)        ! limit thickness to 3.5 m
+      x3=max(x3,0.1)        ! limit thickness >0.1 m since I make fit below
+      x2=min(ICECOEF5*0.5,100.0)  ! convert dia to radius, limit to 100m
+      x2=max(2.5,x2)
+      x2sqr=x2*x2
+      x3sqr=x3*x3
+      amhb = 2.12e-3
+      bmhb = 4.59e-2
+
+      DO IK=1, NK
+        x1=TPI/SIG(IK)   ! period
+        x1sqr=x1*x1
+        KARG1(ik)=-0.26982 + 1.5043*x3 - 0.70112*x3sqr + 0.011037*x2 +  &
+                  (-0.0073178)*x2*x3 + 0.00036604*x2*x3sqr + &
+                  (-0.00045789)*x2sqr + 1.8034e-05*x2sqr*x3 + &
+                  (-0.7246)*x1 + 0.12068*x1*x3 + &
+                  (-0.0051311)*x1*x3sqr + 0.0059241*x1*x2 + &
+                  0.00010771*x1*x2*x3 - 1.0171e-05*x1*x2sqr + &
+                  0.0035412*x1sqr - 0.0031893*x1sqr*x3 + &
+                  (-0.00010791)*x1sqr*x2 + &
+                  0.00031073*x1**3 + 1.5996e-06*x2**3 + 0.090994*x3**3
+        KARG1(IK)=min(KARG1(IK),0.0)
+        ALPHA(IK) = 10.0**KARG1(IK)  
+        perfour=x1sqr*x1sqr
+        if ((x1.gt.5.0) .and. (x1.lt.20.0)) then
+          ALPHA(IK) = ALPHA(IK) + amhb/x1sqr+bmhb/perfour
+        else if (x1.gt.20.0) then
+          ALPHA(IK) = amhb/x1sqr+bmhb/perfour
+        endif
+          WN_I(IK) = ALPHA(IK) * 0.5
+        end do
 
     CASE DEFAULT
       WN_I = ICECOEF1 !Default to IC1: Uniform in k

model/src/w3srcemd.F90

@@ -1,4 +1,4 @@
-!> @file
+
 !> @brief Source term integration routine.
 !>
 !> @author H. L. Tolman
@@ -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
@@ -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
@@ -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
@@ -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
 

model/src/w3str1md.F90

@@ -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 |
@@ -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).
@@ -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)
     !/
@@ -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

model/src/w3wavemd.F90

@@ -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
 
@@ -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.
@@ -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
@@ -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.

regtests/bin/matrix.base

@@ -1957,6 +1957,7 @@
     echo "$rtst                  -w work_IC4_M7                -i input_IC4_M7      $ww3 ww3_tic1.1" >> matrix.body
     echo "$rtst                  -w work_IC4_M8                -i input_IC4_M8      $ww3 ww3_tic1.1" >> matrix.body
     echo "$rtst                  -w work_IC4_M9                -i input_IC4_M9      $ww3 ww3_tic1.1" >> matrix.body
+    echo "$rtst                  -w work_IC4_M10               -i input_IC4_M10     $ww3 ww3_tic1.1" >> matrix.body
     echo "$rtst -g 1000m         -w work_IC5_M1                -i input_IC5_M1      $ww3 ww3_tic1.1" >> matrix.body
     echo "$rtst -g 1000m         -w work_IC5_M2                -i input_IC5_M2      $ww3 ww3_tic1.1" >> matrix.body
     echo "$rtst -g 1000m         -w work_IC5_M3                -i input_IC5_M3      $ww3 ww3_tic1.1" >> matrix.body

regtests/ww3_tic1.1/info

@@ -54,6 +54,7 @@
 #            IC4METHOD = 8 - Meylan et al. (2018) ; Liu et al. (2020)       #
 #                            (NB: redundant with IC5+IC5VEMOD=3)            #
 #            IC4METHOD = 9 - RYW (2021) ; Yu et al. (2022)                  #
+#            IC4METHOD = 10 - Meylan et al. (2021)                          #
 #          IC5 = Choose from three different effective medium models        #
 #            IC5VEMOD  = 1 - Extended Fox and Squire model (EFS)            #
 #            IC5VEMOD  = 2 - Robinson and Palmer model (RP)                 #
@@ -101,6 +102,14 @@
 #               'IC1' 19680606 000000   5.35E-6                             #
 #               'IC2' 19680606 000000   16.05E-6                            #
 #                                                                           #
+#      ------------> &SIC4 IC4METHOD = 10 /                                 #
+#            ...ICECOEF1, ICECOEF5 are required:                            #
+#               T T     Ice parameter 1                                     #
+#               T T     Ice parameter 5                                     #
+#            ...                                                            #
+#               'IC1' 19680606 000000   0.2                                 #
+#               'IC5' 19680606 000000   0.459                               #
+#                                                                           #
 #  Reference (w/plots): Rogers and Orzech, NRL Memorandum Report (2013)     #
 #  available from http://www7320.nrlssc.navy.mil/pubs.php                   #
 #  (This report only covers IC1 and IC2, not IC3, which is newer)           #

regtests/ww3_tic1.1/input_IC4_M10/namelists_1-D.nml

@@ -0,0 +1,2 @@
+&SIC4 IC4METHOD = 10 /
+END OF NAMELISTS

regtests/ww3_tic1.1/input_IC4_M10/points.list

@@ -0,0 +1,16 @@
+0.00 0. 'Point 1 '
+1.00E3 0. 'Point 2 '
+2.00E3 0. 'Point 3 '
+3.00E3 0. 'Point 4 '
+4.00E3 0. 'Point 5 '
+5.00E3 0. 'Point 6 '
+6.00E3 0. 'Point 7 '
+7.00E3 0. 'Point 8 '
+8.00E3 0. 'Point 9 '
+9.00E3 0. 'Point 10 '
+10.00E3 0. 'Point 11 '
+11.00E3 0. 'Point 12 '
+12.00E3 0. 'Point 13 '
+13.00E3 0. 'Point 14 '
+14.00E3 0. 'Point 15 '
+15.00E3 0. 'Point 16 '

regtests/ww3_tic1.1/input_IC4_M10/switch

@@ -0,0 +1 @@
+NOGRB SHRD PR3 UQ FLX2 LN0 ST0 NL0 BT0 DB0 TR0 BS0 IC4 IS0 REF0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7