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

…ST4 code (#1124)

Co-authored-by: Fabrice Ardhuin <[email protected]>

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

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}

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

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, &