correct bugs to run correctly GQM implementation (#1127)

manual/eqs/NL1.tex

@@ -55,37 +55,38 @@ \subsubsection{~$S_{nl}$: Discrete Interaction Approximation (\dia)} \label{sec:
 \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 
+Hence for any $\bk$ one quadruplet selects $\bk_{2,+}$ and $\bk_{3,+}$, and the other quadruplet selects its mirror image 
+$\bk_{2,-}$, $\bk_{2,-}$. Because there are 3 different components interacting in the two DIA-selected quadruplets, any discrete spectral component $(f_r,\theta)$ is actually involved in 6 quadruplets and directly exchanges energy with 12 other components $(f_r',\theta')$. Because the values of $f'_r$ and $\theta'$ do not fall exacly on other discrete components, the spectral density is interpolated using a bilinear interpolation, so that each source term value
+$S_{nl}(\bk)$ contains the direct exchange of energy with 48 other discrete components. 
+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, without making explicit that bilinear interpolation, 
 \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}
+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  
+where the geometry of the quadruplet $(\bk_4,\bk_4,\bk,\bk_5)$ is obtained from that of $(\bk,\bk,\bk_{2,+},\bk_{3,+})$ by a dilation by a factor $(1+\lambda)^2$ and rotation by the angle $\delta_{\theta,2}$;  $(\bk_6,\bk_6,\bk,\bk_7)$ has the same dilation but the opposite rotation; $(\bk_8,\bk_8,\bk_9,\bk)$ is dilated by a factor  $(1-\lambda)^2$ and rotated by the angle $-\delta_{\theta,3}$: and $(\bk_{10},\bk_{10},\bk_{11},\bk)$ is dilated by the same factor and rotated by the opposite angle. 
+
+
+The elementary contributions $\delta  S_{\mathrm{nl}}(\bk_l,\bk_m,\bk_n)$   are given by  
 %----------------------------%
 % Nonlinear interactions DIA %
 %----------------------------%
 % eq:snl_dia
 
 \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] ,
+\delta S_{\mathrm{nl}}(\bk_l,\bk_m,\bk_n) =  \frac{C}{g^4} f_{r,l}^{11} \left [ F_l^2 \left ( \frac{F_m}{(1+\lambda)^4} +
+        \frac{F_n}{(1-\lambda)^4} \right ) - \frac{2 F_l F_m F_n}{(1-\lambda^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.
+where the spectral densities are  $F_l = F(f_{r,l} ,\theta_l)$, 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
 
 \begin{table} \begin{center}
  \begin{tabular}{|l|c|c|} \hline
-                    & $\lambda_{nl}$ &     $C$      \\ \hline
+                    & $\lambda$ &     $C$      \\ \hline
 ST6                 &      0.25      & $3.00 \; 10^7$  \\ \hline
 \wam-3              &      0.25      & $2.78 \; 10^7$  \\ \hline
 ST4 (Ardhuin et al.)&      0.25      & $2.50 \; 10^7$  \\ \hline

manual/manual.bib

@@ -3665,6 +3665,17 @@ @article{art:DC23
 	year = {2023}
 }
 
+@ARTICLE{Webb1978,
+   author = "D. J. Webb",
+   title = "Nonlinear transfer between sea waves",
+   journal = DSR,
+   volume = 25,
+   pages = "279--298",
+   year = 1978,
+   where="paper",
+}
+
+
 @ARTICLE{Lavrenov2001,
    author = "Igor V. Lavrenov",
    title = "Effect of wind wave parameter fluctuation on the nonlinear spectrum evolution",

manual/sys/files_w3.tex

@@ -506,11 +506,15 @@ \subsubsection{~Wave model modules} \label{sec:wave_mod}
 \end{flist}
 
 \noindent
-Nonlinear interaction module (\dia) \hfill {\file w3snl1md.ftn}
+Nonlinear interaction module (\dia or GQM) \hfill {\file w3snl1md.ftn}
 
 \begin{flisti}
 \fit{w3snl1}{Calculation of $S_{nl}$.}
 \fit{insnl1}{Initialization for $S_{nl}$.}
+\fit{w3snlgqm}{Calculation of $S_{nl}$.}
+\fit{w3scouple}{Calculation of coupling coefficient.}
+\fit{gauleg}{Calculation of Gauss-Legendre quadrature coefficients.}
+\fit{INSNLGQM}{Initialization for $S_{nl}$ with GQ method.}
 \end{flisti}
 
 \noindent

model/src/w3snl1md.F90

@@ -825,6 +825,8 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
     USE CONSTANTS, ONLY: TPI
     USE W3GDATMD,  ONLY: SIG, NK ,  NTH , DTH, XFR, FR1, GQTHRSAT, GQAMP
 
+    IMPLICIT NONE
+
     REAL, intent(in) :: A(NTH,NK), CG(NK), WN(NK)
     REAL, intent(in) :: DEPTH
     REAL, intent(out) :: TSTOTn(NTH,NK), TSDERn(NTH,NK)
@@ -883,8 +885,8 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
     !     Gamma_max=1.3 (JFMAX>NF) TO OBTAIN IMPROVED RESULTS
     !     Note by Fabrice Ardhuin: this appears to give the difference in tail benaviour with Gerbrant's WRT
     !=======================================================================
-    JFMIN= 1-INT(LOG(1.0D0)/LOG(RAISF))
-    JFMAX=NF+INT(LOG(1.3D0)/LOG(RAISF))
+    JFMIN=MAX(1-INT(LOG(1.0D0)/LOG(RAISF)),1)
+    JFMAX=MIN(NF+INT(LOG(1.3D0)/LOG(RAISF)),NK)
     !
     !=======================================================================
     !     COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT
@@ -1065,7 +1067,7 @@ SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn)
             TEMP=(TB_TPM(IQ_OM2,JT1,JF1)*(( F(JT1P2P,JFM2)*CF2 *F(JT1P3M,JFM3)*CF3)* &
                  (F(JT,JFM0    )*CF0*TB_V14(JF1)+F(JT1P  ,JFM1)*CF1) &
                  -SP0*SP1P*(SP1P2P*V3_4+SP1P3M*V2_4))+T_2M3P*(AUX05*AUX01-AUX02*AUX06)) *CP0
-            WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT, F(JT,JFM0)
+            WRITE(995,'(3I3,3E12.3)') ICONF,JF,JT, F(JT,JFM0)
             TEMP=(Q_2P3M+Q_2M3P) *CP1
             WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P,  JFM1,AUX00 *CP1, F(JT1P,JFM1),TSTOT(JT1P,JFM1)
             WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P2P,JFM2,-Q_2P3M*CP2,F(JT1P2P,JFM2),TSTOT(JT1P2P,JFM2)
@@ -1219,6 +1221,8 @@ FUNCTION COUPLE(XK1 ,YK1 ,XK2 ,YK2 ,XK3 ,YK3 ,XK4 ,YK4)
     !/ ------------------------------------------------------------------- /
     USE CONSTANTS, ONLY: GRAV
     !
+    IMPLICIT NONE
+
     DOUBLE PRECISION, INTENT(IN)    :: XK1   , YK1   , XK2   , YK2
     DOUBLE PRECISION, INTENT(IN)    :: XK3   , YK3
     DOUBLE PRECISION, INTENT(IN)    :: XK4   , YK4
@@ -1305,6 +1309,7 @@ SUBROUTINE GAULEG (W_LEG ,X_LEG ,NPOIN)
     !/ ------------------------------------------------------------------- /
     !.....VARIABLES IN ARGUMENT
     !     """"""""""""""""""""
+    IMPLICIT NONE
     INTEGER ,         INTENT(IN)    :: NPOIN
     DOUBLE PRECISION ,INTENT(INOUT) :: W_LEG(NPOIN) , X_LEG(NPOIN)
     !
@@ -1552,6 +1557,7 @@ SUBROUTINE INSNLGQM
 #ifdef W3_S
     CALL STRACE (IENT, 'INSNLGQM')
 #endif
+    IMPLICIT NONE
     !.....LOCAL VARIABLES
     INTEGER           JF    , JT    , JF1   , JT1   , NF1P1 , IAUX , NT , NF , IK
     INTEGER           IQ_TE1 , IQ_OM2 , LBUF , DIMBUF , IQ_OM1 , NQ_TE1 , NCONFM
@@ -2084,10 +2090,7 @@ SUBROUTINE INSNLGQM
       AUX=0.0D0
       DO JT1=1,GQNT1
         DO IQ_OM2=1,GQNQ_OM2
-          AAA=TB_FAC(IQ_OM2,JT1,JF1)*TB_TPM(IQ_OM2,JT1,JF1)
-          IF (AAA.GT.AUX) AUX=AAA
-          CCC=TB_FAC(IQ_OM2,JT1,JF1)*TB_TMP(IQ_OM2,JT1,JF1)
-          IF (CCC.GT.AUX) AUX=CCC
+          AUX=MAX(AUX,TB_FAC(IQ_OM2,JT1,JF1)*TB_TPM(IQ_OM2,JT1,JF1),TB_FAC(IQ_OM2,JT1,JF1)*TB_TMP(IQ_OM2,JT1,JF1))
         ENDDO
       ENDDO
       MAXCLA(JF1)=AUX
@@ -2099,6 +2102,7 @@ SUBROUTINE INSNLGQM
     DO JF1=1,GQNF1
       IF (MAXCLA(JF1).GT.AUX) AUX=MAXCLA(JF1)
     ENDDO
+
     TEST1=SEUIL1*AUX
     !
     !.....Set to zero the coupling coefficients not used
@@ -2128,7 +2132,9 @@ SUBROUTINE INSNLGQM
     !
     !..... counts the fraction of the eliminated configurations
     ELIM=(1.D0-DBLE(NCONF)/DBLE(NCONFM))*100.D0
-    !      WRITE(994,*) 'NCONF:',NCONF,ELIM
+#ifdef W3_TGQM
+    WRITE(994,*) 'NCONF, ELIM FRACTION:',NCONF,ELIM
+#endif
   END SUBROUTINE INSNLGQM
   !/
   !/ End of module W3SNL1MD -------------------------------------------- /

model/src/w3src4md.F90

@@ -2520,7 +2520,7 @@ SUBROUTINE W3SDS4 (A, K, CG, USTAR, USDIR, DEPTH, DAIR, SRHS,    &
       RETURN
     END IF
     !
-    WHITECAP(1:2) = 0.
+    WHITECAP(1:4) = 0.
     !
     ! precomputes integration of Lambda over direction
     ! times wavelength times a (a=5 in Reul&Chapron JGR 2003) times dk