remove lookup table for ST4 to speed up computation and clean up the …

…ST4 code (#1124) Co-authored-by: Fabrice Ardhuin <[email protected]>
NOAA-EMC · Jan 16, 2024 · 3952826 · 3952826
1 parent 63f8270
commit 3952826
Show file tree

Hide file tree

Showing 25 changed files with 1,801 additions and 307 deletions.
diff --git a/manual/eqs/ST3.tex b/manual/eqs/ST3.tex
@@ -57,16 +57,15 @@ \subsubsection{~$S_{in} + S_{ds}$: \wam\ cycle 4 (ECWAM)} \label{sec:ST3}
 waves that travel faster than the wind. This accounts for some gustiness in
 the wind and should possibly be resolution-dependent. For reference, this
 parameter was not properly set in early versions of the SWAN model, as
-discovered by R. Lalbeharry.}.  The roughness $z_1$ is defined as,
-
+discovered by R. Lalbeharry.}.  If the friction velocity $u_\star$ is known, 
+it gives the roughness $z_1$ and the wind speed at altitude $z_u$ (by default $z_u=10$~m),  
 \begin{eqnarray}
-U_{10}&=&\frac{u_\star}{\kappa} \log\left(\frac{z_u}{z_1}\right) \\
-z_1&=&\alpha_0 \frac{\tau}{ \sqrt{1-\tau_w/\tau}},
+z_1&=&\alpha_0 \frac{\tau}{ \sqrt{1-\tau_w/u_\star^2}}, \\
+U(z_u)&=&\frac{u_\star}{\kappa} \log\left(\frac{z_u}{z_1}\right) 
 \end{eqnarray}
 
 \noindent
-where $\tau=u_\star^2$, and $z_u$ is the height at which the wind is
-specified. These two equations provide an implicit functional dependence of
+In practice these two equations provide an implicit functional dependence of
 $u_\star$ on $U_{10}$ and $\tau_w/\tau$. This relationship is then tabulated
 \citep{art:Jan91, rep:Bea07}.
 

diff --git a/manual/eqs/ST4.tex b/manual/eqs/ST4.tex
@@ -8,10 +8,10 @@ \subsubsection{~$S_{\mathrm{in}} + S_{\mathrm{ds}}$: Saturation-based dissipatio
 This family of parameterizations uses a positive part of the wind input taken
 from WAM cycle 4 with an ad hoc reduction of $u_\star$, implemented in
 order to allow a balance with a saturation-based dissipation that uses different options for 
-a cumulative term. There are three main options for defining the saturation and the cumulative term. Chosing one or the other is done with the  {\F SDSBCHOICE} parameter, with  {\F SDSBCHOICE=1} for \cite{art:Aea10},  {\F SDSBCHOICE=2} for \cite{Filipot&Ardhuin2012}, and {\F SDSBCHOICE=3} for \cite{Romero2019}. That last options uses a saturation that is defined from the local spectral density, and thus gives zero dissipation for directions where the threshold is not reached, leading to much broader directional spectra. Also the stronger bimodality is achieved by having a strong modulation effect as a cumulative term. 
+a cumulative term. There are three main options for defining the saturation and the cumulative term. Chosing one or the other is done with the  {\F SDSBCHOICE} parameter, with  {\F SDSBCHOICE=1} for \cite{art:Aea10},  {\F SDSBCHOICE=2} for \cite{Filipot&Ardhuin2012}, and {\F SDSBCHOICE=3} for \cite{Romero2019} and later adjustments including \cite{art:AA23}. That last option uses a saturation that is defined from the local spectral density, and thus gives zero dissipation for directions where the threshold is not reached, leading to much broader directional spectra. Also the stronger bimodality is achieved by having a strong modulation effect as a cumulative term. 
 
 Many other adjustments can be made by changing the namelist parameters. A few successful combinations 
-are given by tables \ref{tab:ST4_parSIN} and \ref{tab:ST4_parSDS}, with results described by \citep{art:RA13,art:SAG16}. 
+are given by tables \ref{tab:ST4_parSIN} and \ref{tab:ST4_parSDS}, with results described by \citep{art:RA13,art:SAG16,art:AA23}. 
 Further calibration to any particular wind field should be done for best performance. Guidance for this is given by \cite{Stopa2018}. 
 %We also note that the particular 
 %set of parameters T400 corresponds to setting IPHYS=1 in the ECWAM code cycle 45R2, with a few differences 
@@ -216,27 +216,15 @@ \subsubsection{~$S_{\mathrm{in}} + S_{\mathrm{ds}}$: Saturation-based dissipatio
 direction will typically produce less dissipation than a sea state with all
 the energy radiated in the same direction.
 
-Based on recent analysis by \cite{Guimaraes2018} and \cite{Peureux&al.2019}, this saturation is enhanced by a factor $M_L$ that represents 
-the effect of long waves on short waves 
-\begin{equation}
-M_l(k,\theta)=1+M_\theta \sqrt{\mathrm{mss}(k,\theta)} + N_\theta \sqrt{\mathrm{nss}(k,\theta)} \label{defFACSAT}.
-\end{equation}
-where $M_\theta$ is twice the modulation transfer function for short wave steepness, with 
-$M_\theta=8$ when following the simplified theory by \cite{art:LHS60} and using the root mean square enhancement of $B$ over a 
-long wave cycle. $N_\theta$ is an additional straining factor due to the instability of the wave action envelope of short waves 
-propagating in the direction close to that of the long wave \citep{Peureux&al.2019}. The squared slopes $\mathrm{mss}(k,\theta)$ is 
-the mean square slope in direction $\theta$, wheras $\mathrm{nss}(k,\theta)$ is a slope of long waves propagating in a narrow window $\pm \delta_\theta$, 
-around the short wave direction $\theta$.
-
 We finally define our dissipation term as the sum of the saturation-based term
 and a cumulative breaking term $S_{\mathrm{bk,cu}}$,
 \begin{eqnarray}
 \cS_{ds}(k,\theta)& =&  \sigma
  \frac{C_{\mathrm{ds}}^{\mathrm{sat}}}{B^2_r} \left[ \delta_d
-\max\left\{ M_l(k,\theta) B\left(k\right) -
+\max\left\{  B\left(k\right) -
 B_r,0\right\}^2 \right.
 \nonumber \\
-  & & +  \left(1-\delta_d \right) \left. \max\left\{ M_L(k,\theta) B'\left(k,\theta \right)- B_r
+  & & +  \left(1-\delta_d \right) \left. \max\left\{B'\left(k,\theta \right)- B_r
  ,0\right\}^2\right]N(k,\theta)  \nonumber \\
  & & + \cS_{\mathrm{bk,cu}}(k,\theta) + \cS_{\mathrm{turb}}(k,\theta) \label{Sds_all}.
 \end{eqnarray}

diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90
@@ -897,7 +897,7 @@ MODULE W3GDATMD
     REAL,     POINTER     :: DCKI(:,:), SATWEIGHTS(:,:),CUMULW(:,:),QBI(:,:)
     REAL                  :: AALPHA, BBETA, ZZ0MAX, ZZ0RAT, ZZALP,&
          SSINTHP, TTAUWSHELTER, SSWELLF(1:7), &
-         SSDSC(1:21), SSDSBR,                 &
+         SSDSC(1:21), SSDSBR, SINTAILPAR(1:5),&
          SSDSP, WWNMEANP, SSTXFTF, SSTXFTWN,  &
          FFXPM, FFXFM, FFXFA,   &
          SSDSBRF1, SSDSBRF2, SSDSBINT,SSDSBCK,&
@@ -1317,7 +1317,7 @@ MODULE W3GDATMD
        FFXFM, FFXPM, SSDSBRF1, SSDSBRF2,    &
        SSDSBINT, SSDSBCK, SSDSHCK, SSDSABK, &
        SSDSPBK, SSINBR,SSINTHP,TTAUWSHELTER,&
-       SSWELLF(:), SSDSC(:), SSDSBR,        &
+       SINTAILPAR(:), SSWELLF(:), SSDSC(:), SSDSBR,        &
        SSDSP, WWNMEANP, SSTXFTF, SSTXFTWN,  &
        SSDSBT, SSDSCOS, SSDSDTH, SSDSBM(:)
 #endif
@@ -2074,12 +2074,18 @@ SUBROUTINE W3DIMS  ( IMOD, MK, MTH, NDSE, NDST )
          MPARS(IMOD)%SRCPS%QBI(NKHS,NKD),   &
          STAT=ISTAT                         )
     CHECK_ALLOC_STATUS ( ISTAT )
+    MPARS(IMOD)%SRCPS%IKTAB(:,:)=0.
+    MPARS(IMOD)%SRCPS%DCKI(:,:)=0.
+    MPARS(IMOD)%SRCPS%QBI(:,:)=0.
     SDSNTH  = MTH/2-1 !MIN(NINT(SSDSDTH/(DTH*RADE)),MTH/2-1)
     ALLOCATE( MPARS(IMOD)%SRCPS%SATINDICES(2*SDSNTH+1,MTH), &
          MPARS(IMOD)%SRCPS%SATWEIGHTS(2*SDSNTH+1,MTH), &
          MPARS(IMOD)%SRCPS%CUMULW(MSPEC,MSPEC),        &
          STAT=ISTAT                                   )
     CHECK_ALLOC_STATUS ( ISTAT )
+    MPARS(IMOD)%SRCPS%SATINDICES(:,:)=0.
+    MPARS(IMOD)%SRCPS%SATWEIGHTS(:,:)=0.
+    MPARS(IMOD)%SRCPS%CUMULW(:,:)=0.
 #endif
     !
     SGRDS(IMOD)%SINIT  = .TRUE.
@@ -2648,6 +2654,7 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     ZZ0RAT   => MPARS(IMOD)%SRCPS%ZZ0RAT
     ZZALP    => MPARS(IMOD)%SRCPS%ZZALP
     TTAUWSHELTER  => MPARS(IMOD)%SRCPS%TTAUWSHELTER
+    SINTAILPAR  => MPARS(IMOD)%SRCPS%SINTAILPAR
     SSWELLFPAR  => MPARS(IMOD)%SRCPS%SSWELLFPAR
     SSWELLF  => MPARS(IMOD)%SRCPS%SSWELLF
     SSDSC    => MPARS(IMOD)%SRCPS%SSDSC

diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90
@@ -839,7 +839,8 @@ MODULE W3GRIDMD
 #endif
   !
 #ifdef W3_ST4
-  INTEGER                 :: SWELLFPAR, SDSISO, SDSBRFDF
+  INTEGER                 :: SWELLFPAR, SDSISO, SDSBRFDF, SINTABLE,&
+                             TAUWBUG
   REAL 		   :: SDSBCHOICE
   REAL                    :: ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP,&
        ZALP, Z0RAT, TAUWSHELTER, SWELLF,    &
@@ -855,7 +856,8 @@ MODULE W3GRIDMD
        SDSBRF1,                             &
        SDSBM0, SDSBM1, SDSBM2, SDSBM3,      &
        SDSBM4, SDSFACMTF, SDSCUMP,  SDSNUW, &
-       SDSL, SDSMWD, SDSMWPOW, SPMSS, SDSNMTF
+       SDSL, SDSMWD, SDSMWPOW, SPMSS, SDSNMTF, SINTAIL1, SINTAIL2, &
+       CUMSIGP, VISCSTRESS
 #endif
   !
 #ifdef W3_ST6
@@ -997,7 +999,7 @@ MODULE W3GRIDMD
   NAMELIST /SIN4/ ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP, ZALP, &
        TAUWSHELTER, SWELLFPAR, SWELLF,                 &
        SWELLF2, SWELLF3, SWELLF4, SWELLF5, SWELLF6,    &
-       SWELLF7, Z0RAT, SINBR
+       SWELLF7, Z0RAT, SINBR, SINTABLE, SINTAIL1, SINTAIL2, TAUWBUG, VISCSTRESS
 #endif
 #ifdef W3_NL1
   NAMELIST /SNL1/ LAMBDA, NLPROP, KDCONV, KDMIN,                  &
@@ -1039,7 +1041,7 @@ MODULE W3GRIDMD
        SDSC5, SDSC6, SDSBR, SDSBT, SDSP, SDSISO,       &
        SDSBCK, SDSABK, SDSPBK, SDSBINT, SDSHCK,        &
        SDSDTH, SDSCOS, SDSBRF1, SDSBRFDF,  SDSNUW,     &
-       SDSBM0, SDSBM1, SDSBM2, SDSBM3, SDSBM4,         &
+       SDSBM0, SDSBM1, SDSBM2, SDSBM3, SDSBM4, CUMSIGP,&
        WHITECAPWIDTH, WHITECAPDUR, SDSMWD, SDSMWPOW, SDKOF
 #endif
 
@@ -1718,6 +1720,12 @@ SUBROUTINE W3GRID()
     TAUWSHELTER = 0.3
     ZALP   = 0.006
     SINBR   = 0.
+    SINTABLE = 1
+    SINTAIL1 = 0. !  TAUWSHELTER FOR TAIL (no table)
+    SINTAIL2 = 0. !  additional peak in capillary range
+    TAUWBUG  = 1  !  TAUWBUG is 1 is the bug is kept:
+    !  initializes TAUWX/Y to zero in W3SRCE
+    VISCSTRESS =0
 #endif
     !
 #ifdef W3_ST6
@@ -1801,6 +1809,11 @@ SUBROUTINE W3GRID()
     SSWELLF(6) = SWELLF6
     SSWELLF(7) = SWELLF7
     SSWELLFPAR = SWELLFPAR
+    SINTAILPAR(1) = FLOAT(SINTABLE)
+    SINTAILPAR(2) = SINTAIL1
+    SINTAILPAR(3) = SINTAIL2
+    SINTAILPAR(4) = FLOAT(TAUWBUG)
+    SINTAILPAR(5) = VISCSTRESS
 #endif
     !
 #ifdef W3_ST6
@@ -2106,8 +2119,8 @@ SUBROUTINE W3GRID()
     SDSDTH    = 80.
     SDSCOS    = 2.
     SDSISO    = 2
-    SDSBM0    = 1.
-    SDSBM1    = 0.
+    SDSBM0    = 1.          ! All these parameters are related to finite depth
+    SDSBM1    = 0.          ! scaling of breaking
     SDSBM2    = 0.
     SDSBM3    = 0.
     SDSBM4    = 0.
@@ -2117,8 +2130,9 @@ SUBROUTINE W3GRID()
     SDSBINT   = 0.3
     SDSHCK    = 1.5
     WHITECAPWIDTH = 0.3
-    SDSSTRAIN = 0.
     SDSFACMTF =  400    ! MTF factor for Lambda , Romero (2019)
+    CUMSIGP   = 0.
+    SDSSTRAIN = 0.
     SDSSTRAINA = 15.
     SDSSTRAIN2 = 0.
     WHITECAPDUR   = 0.56 ! breaking duration factor
@@ -2129,7 +2143,7 @@ SUBROUTINE W3GRID()
     ! MTF
     SPMSS     = 0.5    ! cmss^SPMSS
     SDSNMTF   = 1.5    ! MTF power
-    SDSCUMP   = 2.
+    SDSCUMP   = 2.     ! 2 for cumulative mss, 1 for cumulative orb. vel.
     ! MW
     SDSMWD    = .9  ! new AFo
     SDSMWPOW  = 1.  ! (k )^pow
@@ -2211,9 +2225,9 @@ SUBROUTINE W3GRID()
     SSDSC(7)   = WHITECAPWIDTH
     SSDSC(8)   = SDSSTRAIN   ! Straining constant ...
     SSDSC(9)   = SDSL
-    SSDSC(10)  = SDSSTRAINA*NTH/360. ! angle Aor enhanced straining
+    SSDSC(10)  = SDSSTRAINA*NTH/360. ! angle for enhanced straining
     SSDSC(11)  = SDSSTRAIN2  ! straining constant for directional part
-    SSDSC(12)  = SDSBT
+    SSDSC(12)  = CUMSIGP
     SSDSC(13)  = SDSMWD
     SSDSC(14)  = SPMSS
     SSDSC(15)  = SDSMWPOW
@@ -3197,7 +3211,7 @@ SUBROUTINE W3GRID()
 #ifdef W3_ST4
       WRITE (NDSO,2920) ZWND, ALPHA0, Z0MAX, BETAMAX, SINTHP, ZALP,   &
            TAUWSHELTER, SWELLFPAR, SWELLF, SWELLF2, SWELLF3, SWELLF4, &
-           SWELLF5, SWELLF6, SWELLF7, Z0RAT, SINBR
+           SWELLF5, SWELLF6, SWELLF7, Z0RAT, SINBR, SINTABLE, TAUWBUG, VISCSTRESS, SINTAIL1, SINTAIL2
 #endif
 #ifdef W3_ST6
       WRITE (NDSO,2920) SINA0, SINWS, SINFC
@@ -3262,7 +3276,7 @@ SUBROUTINE W3GRID()
            SDSBT, SDSP, SDSISO, SDSCOS, SDSDTH, SDSBRF1,     &
            SDSBRFDF, SDSBM0, SDSBM1, SDSBM2, SDSBM3, SDSBM4, &
            SPMSS, SDKOF, SDSMWD, SDSFACMTF, SDSNMTF,SDSMWPOW,&
-           SDSCUMP, SDSNUW, WHITECAPWIDTH, WHITECAPDUR
+           SDSCUMP, CUMSIGP, SDSNUW, WHITECAPWIDTH, WHITECAPDUR
 #endif
 #ifdef W3_ST6
       WRITE (NDSO,2924) SDSET, SDSA1, SDSA2, SDSP1, SDSP2
@@ -3305,7 +3319,7 @@ SUBROUTINE W3GRID()
            JGS_TERMINATE_DIFFERENCE,                   &
            JGS_TERMINATE_NORM,                         &
            JGS_LIMITER,                                &
-           JGS_LIMITER_FUNC,                           & 
+           JGS_LIMITER_FUNC,                           &
            JGS_USE_JACOBI,                             &
            JGS_BLOCK_GAUSS_SEIDEL,                     &
            JGS_MAXITER,                                &
@@ -3642,7 +3656,7 @@ SUBROUTINE W3GRID()
       END SELECT
 
       IF (FSTOTALIMP .or. FSTOTALEXP) THEN
-        LPDLIB = .TRUE. 
+        LPDLIB = .TRUE.
       ENDIF
       !
       IF (SUM(UNSTSCHEMES).GT.1) WRITE(NDSO,1035)
@@ -6234,7 +6248,9 @@ SUBROUTINE W3GRID()
          '        SWELLF =',F8.5,', SWELLF2 =',F8.5,             &
          ', SWELLF3 =',F8.5,', SWELLF4 =',F9.1,','/              &
          '        SWELLF5 =',F8.5,', SWELLF6 =',F8.5,            &
-         ', SWELLF7 =',F12.2,', Z0RAT =',F8.5,', SINBR =',F8.5,'  /')
+         ', SWELLF7 =',F12.2,', Z0RAT =',F8.5,', SINBR =',F8.5,','/              &
+         '        SINTABLE =',I2,', TAUWBUG =',I2,               &
+         ', VISCSTRESS =',F8.5,', SINTAIL1 =',F8.5,', SINTAIL2 =',F8.5,'  /')
 #endif
     !
 #ifdef W3_ST6
@@ -6407,7 +6423,7 @@ SUBROUTINE W3GRID()
          '        SPMSS = ',F5.2, ', SDKOF =',F5.2,        &
          ', SDSMWD =',F5.2,', SDSFACMTF =',F5.1,', '/      &
          '        SDSMWPOW =',F3.1,', SDSNMTF =', F5.2,    &
-         ', SDSCUMP =', F3.1,', SDSNUW =', E8.3,', '/,     &
+         ', SDSCUMP =', F3.1,', CUMSIGP =', F3.1,', SDSNUW =', E10.3,', '/,     &
          '        WHITECAPWIDTH =',F5.2, ' WHITECAPDUR =',F5.2,' /')
 #endif
     !
@@ -6519,20 +6535,20 @@ SUBROUTINE W3GRID()
 947 FORMAT  (/'  Ice scattering ',A,/ &
          ' --------------------------------------------------')
 948 FORMAT  ('  IS2 Scattering ... '/&
-         '        scattering coefficient       : ',E9.3/ &
-         '        0: no back-scattering        : ',E9.3/ &
+         '        scattering coefficient       : ',E10.3/ &
+         '        0: no back-scattering        : ',E10.3/ &
          '     TRUE: istropic back-scattering  : ',L3/   &
          '     TRUE: update of ICEDMAX         : ',L3/   &
          '     TRUE: keeps updated ICEDMAX     : ',L3/   &
-         '        flexural strength            : ',E9.3/ &
+         '        flexural strength            : ',E10.3/ &
          '     TRUE: uses Robinson-Palmer disp.: ',L3/   &
          '        attenuation                  : ',F5.2/ &
          '        fragility                    : ',F5.2/ &
          '        minimum floe size in meters  : ',F5.2/ &
          '        pack scattering coef 1       : ',F5.2/ &
          '        pack scattering coef 2       : ',F5.2/ &
          '        scaling by concentration     : ',F5.2/ &
-         '        creep B coefficient          : ',E9.3/ &
+         '        creep B coefficient          : ',E10.3/ &
          '        creep C coefficient          : ',F5.2/ &
          '        creep D coefficient          : ',F5.2/ &
          '        creep N power                : ',F5.2/ &
@@ -6543,7 +6559,7 @@ SUBROUTINE W3GRID()
          '        energy of activation         : ',F5.2/ &
          '        anelastic coefficient        : ',E11.3/ &
          '        anelastic exponent           : ',F5.2)
-2948 FORMAT ( '  &SIS2 ISC1 =',E9.3,', IS2BACKSCAT =',E9.3,   &
+2948 FORMAT ( '  &SIS2 ISC1 =',E10.3,', IS2BACKSCAT =',E10.3,   &
          ', IS2ISOSCAT =',L3,', IS2BREAK =',L3,          &
          ', IS2DUPDATE =',L3,','/                        &
          '       IS2FLEXSTR =',E11.3,', IS2DISP =',L3,   &