Merge branch 'develop' into ww3_gr3

manual/defs.tex

@@ -94,6 +94,9 @@
 \newcommand{\cR}{{\cal R}}
 \newcommand{\cS}{{\cal S}}
 
+\newcommand{\rd}{{\mathrm d}}
+
+
 \newcommand{\marbox}[1]{\marginpar{\fbox{{\small #1}}}}
 
 \newcommand{\proddefH}[3]{

manual/eqs/NL1.tex

@@ -1,58 +1,84 @@
-\vsssub
-\subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:NL1}
-\vsssub
-
-\opthead{NL1}{\wam\ model}{H. L. Tolman}
 
 \noindent
-Nonlinear wave-wave interactions can be modeled using the discrete interaction
-approximation \citep[\dia,][]{art:Hea85b}. This parameterization was
+
+
+Resonant nonlinear interactions occur between four wave components
+(quadruplets) with wavenumber vector $\bk$, $\bk_1$, $\bk_2$ and $\bk_3$ are such that 
+% eq:resonance
+\begin{equation} \left .
+\begin{array}{ccc}
+  \bk + \bk_1 & = & \bk_2 + \bk_3     \\
+  f_r + f_{r,1}& =& f_{r,2} +  f_{r,3} 
+\end{array} \:\:\: \right \rbrace \:\:\: , \label{eq:resonance}
+\end{equation}
+
+Nonlinear 4-wave interaction theories were 
 originally developed for the spectrum $F(f_r ,\theta)$. To assure the
 conservative nature of $S_{nl}$ for this spectrum (which can be considered as
 the "final product" of the model), this source term is calculated for
 $F(f_r,\theta)$ instead of $N(k,\theta)$, using the conversion
 (\ref{eq:jac_fr}).
 
-Resonant nonlinear interactions occur between four wave components
-(quadruplets) with wavenumber vector $\bk_1$ through $\bk_4$. In the \dia, it
-is assumed that $\bk_1 = \bk_2$. Resonance conditions then require that
-%--------------------------%
-% Resonance conditions DIA %
-%--------------------------%
+\vsssub
+\subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:NL1}
+\vsssub
+
+\opthead{NL1}{\wam\ model}{H. L. Tolman}
+
+
+
+ In the \dia, for each component $\bk$, only 2 quadruplets configuration are 
+used, while there should be thousands for the full integral, and the interaction caused by these 2 quadruplets 
+is scaled so that it gives the right order of magnitude for the flux of energy towards low frequencies. 
+
+Both quadruplets used the DIA use 
 % eq:resonance
 \begin{equation} \left .
 \begin{array}{ccc}
-  \bk_1 + \bk_2 & = & \bk_3 + \bk_4          \\
-  \sigma_2  & = & \sigma_1                   \\
-  \sigma_3  & = & (1+\lambda_{nl})\sigma_1   \\
-  \sigma_4  & = & (1-\lambda_{nl})\sigma_1
-\end{array} \:\:\: \right \rbrace \:\:\: , \label{eq:resonance}
+  \bk_1 & = & \bk\\
+  f_{r,2}  & = & (1+\lambda)f_{r}    \\
+  f_{r,3}  & = & (1-\lambda)f_{r} 
+\end{array} \:\:\: \right \rbrace \:\:\: , \label{eq:DIAchoice}
 \end{equation}
-where $\lambda_{nl}$ is a constant. For these quadruplets, the contribution
-$\delta S_{nl}$ to the interaction for each discrete $(f_r,\theta)$
-combination of the spectrum corresponding to $\bk_1$ is calculated as
+where $\lambda$ is a constant, usually 0.25, and they only differ by the choice of the interacting angles 
+taking either a plus sign or a minus sign in the following 
+\begin{equation} \left .
+\begin{array}{ccc}
+  \theta_{2,\pm}  & = & \theta \pm \delta_{\theta,2}   \\
+ \theta_{3,\pm}  & = & \theta \mp  \delta_{\theta,3}   \\
+ \end{array} \:\:\: \right \rbrace \:\:\: , \label{eq:DIAangles}
+\end{equation}
+where $\delta_{\theta,2}$ and $\delta_{\theta,3}$ are only a function of $\lambda$ given by the geometry of
+the interacting wavenumbers along the "figure of 8", namely  
+\begin{eqnarray}
+\cos(\delta_{\theta,2})&=&(1-\lambda)^4+4-(1+\lambda)^4)/[4(1-\lambda)^2], \\
+\sin(\delta_{\theta,3})&=&\sin(\delta_{\theta,2}) (1-\lambda)^2/(1+\lambda)^2.
+\end{eqnarray}
+
+ For these quadruplets, each source term value
+$S_{nl}(\bk)$ corresponding to each discrete $(f_r,\theta)$
+we compute the three contributions that correspond to the situation in which $\bk$ takes the role of $\bk$,$\bk_{2,+}$, $\bk_{2,-}$, $\bk_{3,+}$  and $\bk_{3,-}$ in the quadruplet, namely the full source term is 
+\begin{eqnarray} 
+S_{\mathrm{nl}}(\bk) &=& -2  \left[\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,+)+\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,-)\right] \nonumber  \\
+ & & + \delta  S_{\mathrm{nl}}(\bk_4,\bk,\bk_5,+) + \delta  S_{\mathrm{nl}}(\bk_6,\bk,\bk_7,-)  \\
+     & &        +   \delta S_{\mathrm{nl}}(\bk_8,\bk_9,\bk, +) + \delta S_{\mathrm{nl}}(\bk_{10},\bk_{11},\bk, -) . \label{eq:diasum}
+\end{eqnarray}
+with elementary contributions given by  
 %----------------------------%
 % Nonlinear interactions DIA %
 %----------------------------%
 % eq:snl_dia
-\begin{eqnarray}
-\left ( \begin{array}{c}
-  \delta S_{nl,1} \\ \delta S_{nl,3} \\ \delta S_{nl,4}
-\end{array} \right ) & = & D
-\left ( \begin{array}{r} -2 \\ 1 \\ 1 \end{array} \right )
-C g^{-4} f_{r,1}^{11} \times \nonumber \\
-& & \left [ F_1^2
-\left ( \frac{F_3}{(1+\lambda_{nl})^4} +
-        \frac{F_4}{(1-\lambda_{nl})^4} \right ) -
-\frac{2 F_1 F_3 F_4}{(1-\lambda_{nl}^2)^4}
-\right ] \: , \label{eq:snl_dia}
-\end{eqnarray}
-where $F_1 = F(f_{r,1} ,\theta_1 )$ etc. and $\delta S_{nl,1} = \delta
-S_{nl}(f_{r,1} ,\theta_1 )$ etc., $C$ is a proportionality constant. The
-nonlinear interactions are calculated by considering a limited number of
-combinations $(\lambda_{nl},C)$. In practice, only one combination is
-used. Default values for different source term packages are presented in
-Table~\ref{tab:snl_par}.
+
+\begin{equation}
+\delta S_{\mathrm{nl}}(\bk,\bk_2,\bk_3,s) =  \frac{C}{g^4} f_{r,1}^{11} \left [ F^2 \left ( \frac{F_{2,s}}{(1+\lambda_{nl})^4} +
+        \frac{F_{3,s}}{(1-\lambda_{nl})^4} \right ) - \frac{2 F F_{2,s} F_{3,s}}{(1-\lambda_{nl}^2)^4} \right] ,
+      \label{eq:snl_dia}
+\end{equation}
+where $s=+$ or $s=-$ is a sign index, and the spectral densities are  $F = F(f_{r} ,\theta)$, $F_{2,+} = F(f_{r,2} ,\theta + \delta_{\theta,2})$, $F_{2,-} = F(f_{r,2} ,\theta - \delta_{\theta,2})$,  etc.
+ $C$ is a proportionality constant that was tuned to reproduce the inverse energy cascade.  Default values for different source term packages are presented in Table~\ref{tab:snl_par}.
+As a result, when accounting for the two quadruplet configurations, the source term at $\bk$ includes the interactions with 
+10 other spectral components. Besides, because $f_{r,2}$ and $f_{r,3}$ nor $\theta_{2,\pm} $ and $\theta_{3,\pm} $ fall on discretized frequencies and directions, the spectral densities are bilinearly interpolated, which involves 4 discrete spectral components for each of these 10 components. 
+
 
 
 % tab:snl_par
@@ -68,7 +94,7 @@ \subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:
 \caption{Default constants in \dia\ for input-dissipation packages.}
 \label{tab:snl_par} \botline \end{table}
 
-This source term is developed for deep water, using the appropriate dispersion
+This parameterization was developed for deep water, using the appropriate dispersion
 relation in the resonance conditions. For shallow water the expression is
 scaled by the factor $D$ (still using the deep-water dispersion relation,
 however)
@@ -132,3 +158,37 @@ \subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:
 above constants can be reset by the user in the input files of the
 model (see \para\ref{sub:ww3grid}).
 
+\vsssub
+\subsubsection{~$S_{nl}$: Gaussian Quadrature Method  (\dia)} \label{sec:GQM}
+\vsssub
+
+\opthead{NL1 , but with a negative IQTYPE}{TOMAWAC model, M. Benoit}{adaptation to WW3 by S. Siadatmousavi \& F. Ardhuin}
+
+\noindent
+Changing the namelist parameter IQTYPE to a negative value replaces the
+DIA parameterization with the possibility to perform an exact but fast cal-
+culation of $S_{\mathrm{nl}}$ using the Gaussian Quadrature Method of \cite{Lavrenov2001}.
+More details can be found in \cite{Gagnaire-Renou2009}.
+
+
+The quadruplet configurations that are used correspond to the three integrals over $f_1$, $f_2$ and $\theta_1$, with all other frequencies and directions given by the resonance conditions (\ref{eq:resonance}) with only one ambiguity on the angle $\theta_2$ which can be defined by a sign index $s$, as in the DIA.  Starting from eq. (A4) in \cite{Lavrenov2001} as writen in (2.25) of \cite{Gagnaire-Renou2009}, the source term is 
+\begin{equation}
+S_{\mathrm{nl}}(\sigma,\theta) =  8 \sum_s \int_{\sigma_1=0}^\infty \int_{\theta_1=0}^{2 \pi} \int_{\sigma_2=0}^{(\sigma+\sigma_1)/2}  T  \frac{F_2 F_3 (F \sigma_1^4 + F_1 \sigma^4) - F F_1 (F_2 \sigma_3^4 + F_3 \sigma_2^4)}{\sqrt{B}\sqrt{((\left| \bk+\bk_1 \right|/g- \sigma_3^2)^2-\sigma_2^4} } {\mathrm d}\sigma_1 {\mathrm d}\theta_1 {\mathrm d}\sigma_2 ,
+      \label{eq:snl_gqm}
+\end{equation}
+where $B$ is given by eq. (A5) of Lavrenov (2001) and  
+\begin{equation}
+T(\bk,\bk_1,\bk_2,\bk_3) = \frac{\pi g^2 D^2(\bk,\bk_1,\bk_2,\bk_3) }{4 \sigma \sigma_1 \sigma_2 \sigma_3}
+\end{equation}
+where $ D(\bk,\bk_1,\bk_2,\bk_3)$ is given by \cite{Webb1978} in his eq. (A1). 
+
+This triple integral is performed using quadrature functions to best resolve the effect of the singularities in the denominator. It is thus replaced with weighted sums over the 3 dimensions. 
+
+Compared to the DIA, there is no bilinear interpolation and the nearest neighbor is used in frequency and direction. Also, 
+the source term is computed by a loop over the quadruplet configuration, which allows for filtering based on 
+both the value of the coupling coefficient and the energy level at the frequency corresponding to $\bk$. Within 
+that loop, the source term contribution is computed for all 4 interacting components, so that any filtering still 
+conserves energy, action, momentum ... (One may argue that this multiplies by 4 the number of calculations, but it may have the benefit of properly dealing with the high frequency boundary... this is to be verified. The same question arises for the DIA: why have the wavenumber $\bk$ play the role of the other members of the quadruplets when this will also be computed as we loop on the spectral components?). 
+
+If a very aggressive filtering is performed, the source may need to be rescaled.
+

manual/eqs/output.tex

@@ -12,9 +12,9 @@ \subsection{~Output parameters} \label{sub:outpars}
 in \para\ref{sec:ww3shel}. That input file also provides a list of flags  
 indicating if output parameters are available in different field
 output file types (ASCII, grib, igrads, NetCDF). 
-For any details on how these parameters are computed, the user may read the code of the  {\code w3iogo} routine, in the {\code w3iogomd.ftn} module. 
+For any details on how these parameters are computed, the user may read the code of the  {\code w3iogo} routine, in the {\code w3iogomd.F90} module. 
 
-Selection of field outputs in  {\code ww3\_shel.inp} is most easily performed by providing a list of the 
+Selection of field outputs in {\code ww3\_shel.nml} or {\code ww3\_shel.inp} is most easily performed by providing a list of the 
 requested parameters, for example, {\textbf HS DIR SPR}  will request the calculation of significant wave height, mean direction and directional spread. These will thus be stored in the {\code out\_grd.XX} file and can be post-processed, for example in NetCDF using   {\code ww3\_ouf}. Examples are given in \para\ref{sec:ww3multi} and
 \para\ref{sec:ww3ounf}. The names for these namelists are the bold names below, for 
 example \textbf{HS}. 
@@ -26,6 +26,9 @@ \subsection{~Output parameters} \label{sub:outpars}
 file extensions, NetCDF variable names and namelist-based selection (see 
 also \para\ref{sec:ww3ounf}), and the long parameter name/definition.
 
+When the result is not overly sensitive to the contribution of the unresolved part of the spectrum (for $f<f_{NK}$), the contribution of the tail is parameterized assuming a power law decay of the spectrum, by default $E(f,\theta) = E(f_{NK},\theta) (f_{NK}/f)^{-5}$, for some parameters this is either unnecessary or misleading, and the integrals are computed only up to $f_{NK}$. 
+
+
 Finally we note that in all definitions the frequency is the \emph{relative} frequency. Thus, in the presence of currents, these can only be compared to drifting measurement data or data obtained in wavenumber and converted to frequency. Comparison to fixed instrument data requires the use of the full spectrum and proper conversion to the fixed reference frame.
 
 \begin{list}{\Roman{outgrps})\hfill}
@@ -300,8 +303,12 @@ \subsection{~Output parameters} \label{sub:outpars}
 \item \textbf{MSD} Direction of the maximum slope variance mss$_u$
 \item \textbf{MCD} Spectral tail direction
 \item \textbf{QP}  Peakedness parameter \citep{art:G70}
-      \begin{equation} Q_p = \frac{2}{E^2} \int_0^{2\pi} \int_0^\infty
-      \sigma\:F(\sigma,\theta)^2\:d\sigma\:d\theta \: \label{eq:qp}
+      \begin{equation} Q_p = \frac{2}{E^2} \int_0^{f_{NK}} f \left( \int_0^{2\pi} 
+      F(f,\theta) \:\rd \theta \right)^2 \: \rd  f \: \label{eq:qp}
+      \end{equation}
+\item \textbf{QKK}  wavenumber peakedness \citep{art:DC23}
+      \begin{equation} Q_{kk} = \frac{1}{E^2} \int_0^{f_{NK}} \int_0^{2\pi} 
+      0.5 \left[ A(k,\theta)+ A(k,\theta+\pi)\right]^2 \frac{\sigma^2}{k C_g} \:\rd \theta  \: \rd  \sigma \: \label{eq:qkk}
       \end{equation}
 \end{list}
 

manual/impl/switch.tex

@@ -94,7 +94,7 @@ \subsubsection{~Mandatory switches} \label{sub:man_switch}
 Selection of nonlinear interactions:
 \begin{slist}
 \sit{nl0} {No nonlinear interactions used.}
-\sit{nl1} {Discrete interaction approximation (\dia).}
+\sit{nl1} {Discrete interaction approximation (\dia) or Gaussian Quadrature Method (GQM).}
 \sit{nl2} {Exact interaction approximation (\xnl).}
 \sit{nl3} {Generalized Multiple \dia\ (\gmd).}
 \sit{nl4} {Two-scale approximation (TSA).} 

manual/manual.bib

@@ -3654,3 +3654,33 @@ @article{vanVledder2006
 	volume = {53},
 	year = {2006}
 }
+
+@article{art:DC23,
+	author = {De Carlo, Marine and Fabrice Ardhuin and Annabelle Ollivier and Adrien Nigou },
+	journal = {Journal of Geophysical Research - Oceans},
+	keywords = {wave groups,altimetry},
+	pages = {},
+	title = {Wave groups and small scale variability of wave heights observed by altimeters},
+	volume = {},
+	year = {2023}
+}
+
+@ARTICLE{Lavrenov2001,
+   author = "Igor V. Lavrenov",
+   title = "Effect of wind wave parameter fluctuation on the nonlinear spectrum evolution",
+   journal = JPO,
+   volume = 31,
+   pages = "861--873",
+   year = 2001,
+   url="http://ams.allenpress.com/archive/1520-0485/31/4/pdf/i1520-0485-31-4-861",
+   keywords={4-wave interactions,GQM},
+}
+
+
+@PHDTHESIS{Gagnaire-Renou2009,
+   author = "Elodie Gagnaire-Renou",
+   title = "Amelioration de la modelisation spectrale des etats de mer par un calcul quasi-exact des interactions non-lineaires vague-vague",
+   school = "Universit{\'e} du Sud Toulon Var",
+   year = 2010,
+}
+

model/inp/ww3_grid.inp

@@ -102,6 +102,18 @@ $                           KDCONV : Factor before kd in Eq. (n.nn).
 $                           KDMIN, SNLCS1, SNLCS2, SNLCS3 :
 $                                    Minimum kd, and constants c1-3
 $                                    in depth scaling function.
+$                           IQTYPE : Type of depth treatment
+$                                    -2 : Deep water GQM with scaling
+$                                     1 : Deep water DIA
+$                                     2 : Deep water DIA with scaling
+$                                     3 : Shallow water DIA
+$                           TAILNL : Parametric tail power.
+$                           GQMNF1 : number of points along the locus
+$                           GQMNT1 : idem
+$                           GQMNQ_OM2 : idem
+$                           GQMTHRSAT : threshold on saturation for SNL calculation
+$                           GQMTHRCOU : threshold for filter on coupling coefficient
+$                           GQAMP1, GQAMP2, GQAMP3, GQAMP4 : amplification factor
 $   Exact interactions  : Namelist SNL2
 $                           IQTYPE : Type of depth treatment
 $                                     1 : Deep water

model/inp/ww3_shel.inp

@@ -139,6 +139,7 @@ $  T  T  2    15   HMAXD      SDMH  St Dev of MXC (STE)
 $  T  T  2    16   HCMAXD     SDMHC St Dev of MXHC (STE)
 $  F  T  2    17   WBT        WBT   Dominant wave breaking probability bT
 $  F  F  2    18   FP0        TP    Peak period (from peak freq)
+$  F  F  2    19   WNMEAN     WNM   Mean wavenumber
 $   -------------------------------------------------
 $        3                          Spectral Parameters (first 5)
 $   -------------------------------------------------
@@ -211,10 +212,10 @@ $        8                          Spectrum parameters
 $   -------------------------------------------------
 $  F  F  8     1   MSS[X,Y]   MSS   Mean square slopes
 $  F  F  8     2   MSC[X,Y]   MSC   Spectral level at high frequency tail
-$  F  F  8     3   WL02[X,Y]  WL02  East/X North/Y mean wavelength compon
-$  F  F  8     4   ALPXT      AXT   Correl sea surface gradients (x,t)
-$  F  F  8     5   ALPYT      AYT   Correl sea surface gradients (y,t)
-$  F  F  8     6   ALPXY      AXY   Correl sea surface gradients (x,y)
+!  F  F  8     3   MSSD       MSD   Slope direction
+!  F  F  8     4   MSCD       MCD   Tail slope direction
+!  F  F  8     5   QP         QP    Goda peakedness parameter
+!  F  F  8     6   QKK        QKK   Wavenumber peakedness
 $   -------------------------------------------------
 $        9                          Numerical diagnostics
 $   -------------------------------------------------

model/nml/namelists.nml

@@ -81,6 +81,16 @@ $                           KDCONV : Factor before kd in Eq. (n.nn).
 $                           KDMIN, SNLCS1, SNLCS2, SNLCS3 :
 $                                    Minimum kd, and constants c1-3
 $                                    in depth scaling function.
+$                           IQTYPE : Type of depth treatment
+$                                    -2 : Deep water GQM with scaling
+$                                     1 : Deep water DIA
+$                                     2 : Deep water DIA with scaling
+$                                     3 : Shallow water DIA
+$                           TAILNL : Parametric tail power.
+$                           GQMNF1, GQMNT1, GQMNQ_OM2 : Gaussian quadrature resolution
+$                           GQMTHRSAT : Threshold on saturation for SNL calculation
+$                           GQMTHRCOU : Threshold for filter on coupling coefficient
+$                           GQAMP1, GQAMP2, GQAMP3, GQAMP4 : Amplification factors
 $   Exact interactions  : Namelist SNL2
 $                           IQTYPE : Type of depth treatment
 $                                     1 : Deep water

model/nml/ww3_multi.nml

@@ -195,13 +195,13 @@
 !
 ! * the detailed list of field names is given in model/nml/ww3_shel.nml :
 !  DPT CUR WND AST WLV ICE IBG TAU RHO D50 IC1 IC5
-!  HS LM T02 T0M1 T01 FP DIR SPR DP HIG
+!  HS LM T02 T0M1 T01 FP DIR SPR DP HIG MXE MXES MXH MXHC SDMH SDMHC WBT TP WNM
 !  EF TH1M STH1M TH2M STH2M WN
 !  PHS PTP PLP PDIR PSPR PWS PDP PQP PPE PGW PSW PTM10 PT01 PT02 PEP TWS PNR
 !  UST CHA CGE FAW TAW TWA WCC WCF WCH WCM FWS
 !  SXY TWO BHD FOC TUS USS P2S USF P2L TWI FIC USP TOC
 !  ABR UBR BED FBB TBB
-!  MSS MSC WL02 AXT AYT AXY
+!  MSS MSC MSD MCD QP QKK
 !  DTD FC CFX CFD CFK
 !  U1 U2
 !