From 8ff8d22f1d76b4cb60a62508e5637ecda4b7188c Mon Sep 17 00:00:00 2001 From: ma05cde Date: Tue, 21 Nov 2023 14:11:38 +0000 Subject: [PATCH 1/2] correct bugs for GQM to run correctly --- manual/eqs/NL1.tex | 29 +++++++++++++++-------------- manual/manual.bib | 11 +++++++++++ manual/sys/files_w3.tex | 6 +++++- model/src/w3snl1md.F90 | 22 ++++++++++++++-------- 4 files changed, 45 insertions(+), 23 deletions(-) diff --git a/manual/eqs/NL1.tex b/manual/eqs/NL1.tex index d9bc4c521..c45e9294d 100644 --- a/manual/eqs/NL1.tex +++ b/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 diff --git a/manual/manual.bib b/manual/manual.bib index 0aacb5105..c49e3340b 100644 --- a/manual/manual.bib +++ b/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", diff --git a/manual/sys/files_w3.tex b/manual/sys/files_w3.tex index fcd48a8f7..d0ad76f7e 100644 --- a/manual/sys/files_w3.tex +++ b/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 diff --git a/model/src/w3snl1md.F90 b/model/src/w3snl1md.F90 index 598b627ea..09c096d2b 100644 --- a/model/src/w3snl1md.F90 +++ b/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 -------------------------------------------- / From 77619fe294054dea62d78f277c480e3b6095894a Mon Sep 17 00:00:00 2001 From: ma05cde Date: Tue, 21 Nov 2023 14:21:14 +0000 Subject: [PATCH 2/2] bugfix on WHITECAP initialization --- model/src/w3src4md.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/model/src/w3src4md.F90 b/model/src/w3src4md.F90 index e2bf12c9a..a1d4423bf 100644 --- a/model/src/w3src4md.F90 +++ b/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