implementation of the GQM (Gaussian Quadrature Method) to replace the…

… DIA in NL1 or NL2. (NOAA-EMC#1083)

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/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

@@ -3664,3 +3664,23 @@ @article{art:DC23
 	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/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/src/w3gdatmd.F90

@@ -429,6 +429,17 @@ MODULE W3GDATMD
   !      KDCON     Real  Public   Conversion factor for relative depth.
   !      KDMN      Real  Public   Minimum relative depth.
   !      SNLSn     Real  Public   Constants in shallow water factor.
+  !      IQTPE     Int.  Public   Type of depth treatment
+  !                               -2 : Deep water GQM with scaling
+  !                                1 : Deep water DIA
+  !                                2 : Deep water DIA with scaling
+  !                                3 : Finite water depth DIA
+  !      GQNF1     Int.  Public   Gaussian quadrature resolution
+  !      GQNT1     Int.  Public   Gaussian quadrature resolution
+  !      GQNNQ_OM2 Int.  Public   Gaussian quadrature resolution
+  !      GQTHRSAT  Real  Public   Threshold on saturation for SNL calculation
+  !      GQTHRCOU  Real  Public   Threshold for filter on coupling coefficient
+  !      GQAMP     R.A.  Public   Amplification factors
   !                                                             (!/NL2)
   !      IQTPE     Int.  Public   Type of depth treatment
   !                                1 : Deep water
@@ -910,6 +921,8 @@ MODULE W3GDATMD
 #ifdef W3_NL1
     REAL                  :: SNLC1, LAM, KDCON, KDMN,             &
          SNLS1, SNLS2, SNLS3
+    INTEGER               :: IQTPE, GQNF1, GQNT1, GQNQ_OM2
+    REAL                  :: NLTAIL, GQTHRSAT, GQTHRCOU, GQAMP(4)
 #endif
 #ifdef W3_NL2
     INTEGER               :: IQTPE, NDPTHS
@@ -1319,6 +1332,8 @@ MODULE W3GDATMD
   !/ Data aliasses for structure SNLP(S)
   !/
 #ifdef W3_NL1
+  INTEGER, POINTER        :: IQTPE, GQNF1, GQNT1, GQNQ_OM2
+  REAL, POINTER           :: NLTAIL, GQTHRSAT, GQTHRCOU, GQAMP(:)
   REAL, POINTER           :: SNLC1, LAM, KDCON, KDMN,             &
        SNLS1, SNLS2, SNLS3
 #endif
@@ -2690,6 +2705,14 @@ SUBROUTINE W3SETG ( IMOD, NDSE, NDST )
     SNLS1  => MPARS(IMOD)%SNLPS%SNLS1
     SNLS2  => MPARS(IMOD)%SNLPS%SNLS2
     SNLS3  => MPARS(IMOD)%SNLPS%SNLS3
+    IQTPE  => MPARS(IMOD)%SNLPS%IQTPE
+    GQNF1  => MPARS(IMOD)%SNLPS%GQNF1
+    GQNT1  => MPARS(IMOD)%SNLPS%GQNT1
+    GQNQ_OM2  => MPARS(IMOD)%SNLPS%GQNQ_OM2
+    NLTAIL => MPARS(IMOD)%SNLPS%NLTAIL
+    GQTHRSAT => MPARS(IMOD)%SNLPS%GQTHRSAT
+    GQTHRCOU=> MPARS(IMOD)%SNLPS%GQTHRCOU
+    GQAMP=> MPARS(IMOD)%SNLPS%GQAMP
 #endif
 #ifdef W3_NL2
     IQTPE  => MPARS(IMOD)%SNLPS%IQTPE