diff --git a/manual/defs.tex b/manual/defs.tex index 1f67da7e7..5b8963fa5 100644 --- a/manual/defs.tex +++ b/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]{ diff --git a/manual/eqs/NL1.tex b/manual/eqs/NL1.tex index a6539dbe5..d9bc4c521 100644 --- a/manual/eqs/NL1.tex +++ b/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. + diff --git a/manual/impl/switch.tex b/manual/impl/switch.tex index 856a7d218..22ab75f34 100644 --- a/manual/impl/switch.tex +++ b/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).} diff --git a/manual/manual.bib b/manual/manual.bib index 8cb9734fb..0aacb5105 100644 --- a/manual/manual.bib +++ b/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, +} + diff --git a/model/inp/ww3_grid.inp b/model/inp/ww3_grid.inp index b802c6717..655a10493 100644 --- a/model/inp/ww3_grid.inp +++ b/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 diff --git a/model/nml/namelists.nml b/model/nml/namelists.nml index 7b373c71e..390fdb874 100644 --- a/model/nml/namelists.nml +++ b/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 diff --git a/model/src/w3gdatmd.F90 b/model/src/w3gdatmd.F90 index 6cd6e91d8..7bc5e2f30 100644 --- a/model/src/w3gdatmd.F90 +++ b/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 diff --git a/model/src/w3gridmd.F90 b/model/src/w3gridmd.F90 index aa618b59f..fa8128afb 100644 --- a/model/src/w3gridmd.F90 +++ b/model/src/w3gridmd.F90 @@ -113,6 +113,7 @@ MODULE W3GRIDMD !/ 27-May-2021 : Moved to a subroutine ( version 7.13 ) !/ 07-Jun-2021 : S_{nl} GKE NL5 (Q. Liu) ( version 7.13 ) !/ 19-Jul-2021 : Momentum and air density support ( version 7.14 ) + !/ 28-Feb-2023 : GQM as an alternative for NL1 ( version 7.15 ) !/ !/ Copyright 2009-2013 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -439,7 +440,7 @@ MODULE W3GRIDMD ! (2006) input and Babanin et al. (2001,2010) dissipation. ! ! !/NL0 No nonlinear interactions. - ! !/NL1 Discrete interaction approximation (DIA). + ! !/NL1 Discrete interaction approximation (DIA or GQM). ! !/NL2 Exact interactions (WRT). ! !/NL3 Generalized Multiple DIA (GMD). ! !/NL4 Two Scale Approximation @@ -867,6 +868,8 @@ MODULE W3GRIDMD #ifdef W3_NL1 REAL :: LAMBDA, KDCONV, KDMIN, & SNLCS1, SNLCS2, SNLCS3 + INTEGER :: IQTYPE, GQMNF1, GQMNT1, GQMNQ_OM2 + REAL :: TAILNL, GQMTHRSAT, GQMTHRCOU, GQAMP1, GQAMP2, GQAMP3, GQAMP4 #endif #ifdef W3_NL2 INTEGER :: IQTYPE, NDEPTH @@ -998,7 +1001,9 @@ MODULE W3GRIDMD #endif #ifdef W3_NL1 NAMELIST /SNL1/ LAMBDA, NLPROP, KDCONV, KDMIN, & - SNLCS1, SNLCS2, SNLCS3 + SNLCS1, SNLCS2, SNLCS3, & + IQTYPE, TAILNL, GQMNF1, GQMNT1, & + GQMNQ_OM2, GQMTHRSAT, GQMTHRCOU, GQAMP1, GQAMP2, GQAMP3, GQAMP4 #endif #ifdef W3_NL2 NAMELIST /SNL2/ IQTYPE, TAILNL, NDEPTH @@ -1831,6 +1836,18 @@ SUBROUTINE W3GRID() SNLCS1 = 5.5 SNLCS2 = 0.833 SNLCS3 = -1.25 + ! Additional parameters for GQM + IQTYPE = 1 + TAILNL = -FACHF + GQMNF1 = 14 + GQMNT1 = 8 + GQMNQ_OM2=8 + GQMTHRSAT=0. + GQMTHRCOU=0.015 + GQAMP1=1. + GQAMP2=0.002 + GQAMP3=1. + GQAMP4=1. CALL READNL ( NDSS, 'SNL1', STATUS ) WRITE (NDSO,922) STATUS WRITE (NDSO,923) LAMBDA, NLPROP, KDCONV, KDMIN, & @@ -1842,6 +1859,18 @@ SUBROUTINE W3GRID() SNLS1 = SNLCS1 SNLS2 = SNLCS2 SNLS3 = SNLCS3 + ! Additional parameters for GQM + IQTPE = IQTYPE + GQNF1 = GQMNF1 + GQNT1 = GQMNT1 + GQNQ_OM2 = GQMNQ_OM2 + GQTHRSAT = GQMTHRSAT + GQTHRCOU = GQMTHRCOU + GQAMP(1) = GQAMP1 + GQAMP(2) = GQAMP2 + GQAMP(3) = GQAMP3 + GQAMP(4) = GQAMP4 + NLTAIL = TAILNL #endif ! #ifdef W3_ST0 @@ -3175,7 +3204,10 @@ SUBROUTINE W3GRID() #endif #ifdef W3_NL1 WRITE (NDSO,2922) LAMBDA, NLPROP, KDCONV, KDMIN, & - SNLCS1, SNLCS2, SNLCS3 + SNLCS1, SNLCS2, SNLCS3, & + IQTYPE, TAILNL, GQMNF1, & + GQMNT1, GQMNQ_OM2, GQMTHRSAT, GQMTHRCOU,& + GQAMP1, GQAMP2, GQAMP3, GQAMP4 #endif #ifdef W3_NL2 WRITE (NDSO,2922) IQTYPE, TAILNL, NDEPTH @@ -6230,7 +6262,11 @@ SUBROUTINE W3GRID() 2922 FORMAT ( ' &SNL1 LAMBDA =',F7.3,', NLPROP =',E10.3, & ', KDCONV =',F7.3,', KDMIN =',F7.3,','/ & ' SNLCS1 =',F7.3,', SNLCS2 =',F7.3, & - ', SNLCS3 = ',F7.3,' /') + ', SNLCS3 = ',F7.3','/ & + ' IQTYPE =',I2,', TAILNL =',F5.1,','/ & + ' GQMNF1 =',I2,', GQMNT1 =',I2,',', & + ' GQMNQ_OM2 =',I2,', GQMTHRSAT =',E11.4,', GQMTHRCOU =',F4.3,','/ & + ' GQAMP1 =',F5.3,', GQAMP2 =',F5.3,', GQAMP3 =',F5.3,', GQAMP4 =',F5.3,' /') #endif ! #ifdef W3_NL2 diff --git a/model/src/w3iogrmd.F90 b/model/src/w3iogrmd.F90 index e96b7b788..4f211402d 100644 --- a/model/src/w3iogrmd.F90 +++ b/model/src/w3iogrmd.F90 @@ -284,7 +284,7 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT & DIKCUMUL #endif #ifdef W3_NL1 - USE W3SNL1MD, ONLY: INSNL1 + USE W3SNL1MD, ONLY: INSNL1, INSNLGQM #endif #ifdef W3_NL2 USE W3SNL2MD, ONLY: INSNL2 @@ -1580,18 +1580,28 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT & #ifdef W3_NL1 IF ( WRITE ) THEN WRITE (NDSM) & - SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3 + SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3, & + IQTPE, NLTAIL, GQNF1, GQNT1, & + GQNQ_OM2, GQTHRSAT, GQTHRCOU, GQAMP #ifdef W3_ASCII WRITE (NDSA,*) & - 'SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3:',& - SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3 + 'SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3, & + IQTPE, NLTAIL, GQNF1, GQNT1, & + GQNQ_OM2, GQTHRSAT, GQTHRCOU, GQAMP:', & + SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3, & + IQTPE, NLTAIL, GQNF1, GQNT1, & + GQNQ_OM2, GQTHRSAT, GQTHRCOU, GQAMP #endif ELSE READ (NDSM,END=801,ERR=802,IOSTAT=IERR) & - SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3 + SNLC1, LAM, KDCON, KDMN, SNLS1, SNLS2, SNLS3, & + IQTPE, NLTAIL, GQNF1, GQNT1, & + GQNQ_OM2, GQTHRSAT, GQTHRCOU, GQAMP END IF IF ( FLTEST ) WRITE (NDST,9051) SNLC1, LAM, & - KDCON, KDMN, SNLS1, SNLS2, SNLS3 + KDCON, KDMN, SNLS1, SNLS2, SNLS3, & + IQTPE, NLTAIL, GQNF1, GQNT1, GQNQ_OM2, & + GQTHRSAT, GQTHRCOU, GQAMP #endif ! #ifdef W3_NL2 @@ -1713,7 +1723,13 @@ SUBROUTINE W3IOGR ( INXOUT, NDSM, IMOD, FEXT & #endif ! #ifdef W3_NL1 - IF ( .NOT. WRITE ) CALL INSNL1 ( IGRD ) + IF ( .NOT. WRITE ) THEN + IF (IQTPE.GT.0) THEN + CALL INSNL1 ( IGRD ) + ELSE + CALL INSNLGQM + END IF + END IF #endif #ifdef W3_NL3 IF ( .NOT. WRITE ) CALL INSNL3 diff --git a/model/src/w3snl1md.F90 b/model/src/w3snl1md.F90 index e21349ede..598b627ea 100644 --- a/model/src/w3snl1md.F90 +++ b/model/src/w3snl1md.F90 @@ -28,7 +28,7 @@ MODULE W3SNL1MD !/ | WAVEWATCH III NOAA/NCEP | !/ | H. L. Tolman | !/ | FORTRAN 90 | - !/ | Last update : 03-Sep-2012 | + !/ | Last update : 28-Feb-2023 | !/ +-----------------------------------+ !/ !/ 04-Feb-2000 : Origination. ( version 2.00 ) @@ -36,6 +36,7 @@ MODULE W3SNL1MD !/ 24-Dec-2004 : Multiple grid version. ( version 3.06 ) !/ 29-May-2009 : Preparing distribution version. ( version 3.14 ) !/ 03-Sep-2012 : Clean up of test output T0, T1 ( version 4.07 ) + !/ 28-Feb-2023 : Adds GQM separate routines ( version 7.07 ) !/ !/ Copyright 2009 National Weather Service (NWS), !/ National Oceanic and Atmospheric Administration. All rights @@ -77,8 +78,22 @@ MODULE W3SNL1MD ! !/ ------------------------------------------------------------------- / !/ + !/ PUBLIC !/ + !/ These are the arrays and variables used for GQM method + !/ + INTEGER :: NCONF + INTEGER, ALLOCATABLE :: K_IF2 (:,:,:) , K_IF3 (:,:,:) , K_1P2P(:,:,:) , & + K_1P3M(:,:,:) , K_1P2M(:,:,:) , K_1P3P(:,:,:) , & + K_1M2P(:,:,:) , K_1M3M(:,:,:) , K_1M2M(:,:,:) , & + K_1M3P(:,:,:) + INTEGER, ALLOCATABLE :: F_POIN(:) , T_POIN(:) , K_IF1(:) , K_1P(:,:) , & + K_1M(:,:) , IDCONF(:,:) + DOUBLE PRECISION, ALLOCATABLE :: F_COEF(:) , F_PROJ(:) , TB_SCA(:) , TB_V14(:) + DOUBLE PRECISION, ALLOCATABLE :: TB_V24(:,:,:) , TB_V34(:,:,:) , & + TB_TPM(:,:,:) , TB_TMP(:,:,:) , TB_FAC(:,:,:) + !/ CONTAINS !/ ------------------------------------------------------------------- / @@ -768,6 +783,1353 @@ SUBROUTINE INSNL1 ( IMOD ) !/ End of INSNL1 ----------------------------------------------------- / !/ END SUBROUTINE INSNL1 + + !/ ------------------------------------------------------------------- / + SUBROUTINE W3SNLGQM(A,CG,WN,DEPTH,TSTOTn,TSDERn) + ! This and the following routines are adapted to WW3 from TOMAWAC qnlin3.f + !*********************************************************************** + ! TOMAWAC V6P1 24/06/2011 + !*********************************************************************** + ! + !brief COMPUTES THE CONTRIBUTION OF THE NON-LINEAR INTERACTIONS + !+ SOURCE TERM BETWEEN QUADRUPLETS USING THE GQM METHOD + !+ ("GAUSSIAN QUADRATURE METHOD") PROPOSED BY LAVRENOV + !+ (2001) + !+ + !+ PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES + !+ FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS + !+ ARE EVENLY DISTRIBUTED OVER [0;2.PI]. + ! + !note THIS SUBROUTINE USES THE OUTPUT FROM 'PRENL3' TO OPTIMISE + !+ THE COMPUTATIONS FOR DIA. + ! + !reference LAVRENOV, I.V. (2001): + !+ "EFFECT OF WIND WAVE PARAMETER FLUCTUATION ON THE NONLINEAR + !+ SPECTRUM EVOLUTION". J. PHYS. OCEANOGR. 31, 861-873. + ! + !history E. GAGNAIRE-RENOU + !+ 04/2011 + !+ V6P1 + !+ CREATED + ! + !history G.MATTAROLO (EDF - LNHE) + !+ 24/06/2011 + !+ V6P1 + !+ Translation of French names of the variables in argument + + ! + !/ Warning, contrary to the DIA routine, there is no extension to frequencies below IK=1 + !/ as a result the first two frequencies are not fully treated. + !================================================================================== + ! This subroutine is same as qnlin3 in TOMWAC + USE CONSTANTS, ONLY: TPI + USE W3GDATMD, ONLY: SIG, NK , NTH , DTH, XFR, FR1, GQTHRSAT, GQAMP + + REAL, intent(in) :: A(NTH,NK), CG(NK), WN(NK) + REAL, intent(in) :: DEPTH + REAL, intent(out) :: TSTOTn(NTH,NK), TSDERn(NTH,NK) + + INTEGER :: ITH,IK,NT,NF + REAL :: q_dfac, SATVAL(NK), SUME, ACCVAL, ACCMAX, AMPFAC + DOUBLE PRECISION :: RAISF, FREQ(NK) + DOUBLE PRECISION :: TSTOT(NTH,NK) , TSDER(NTH,NK), F(NTH,NK) + DOUBLE PRECISION :: TEMP + + !.....LOCAL VARIABLES + INTEGER JF , JT , JF1 , JT1 , IQ_OM2 & + , JFM0 , JFM1 , JFM2 , JFM3 , IXF1 , IXF2 & + , IXF3 , JFMIN , JFMAX , ICONF , LBUF + INTEGER KT1P , KT1M , JT1P , JT1M , KT1P2P, KT1P2M & + , KT1P3P, KT1P3M, KT1M2P, KT1M2M, KT1M3P, KT1M3M & + , JT1P2P, JT1P2M, JT1P3P, JT1P3M, JT1M2P, JT1M2M & + , JT1M3P, JT1M3M + DOUBLE PRECISION V1_4 , V2_4 , V3_4 , Q_2P3M, Q_2M3P, FACTOR & + , T_2P3M, T_2M3P, S_2P3M, S_2M3P, SCAL_T, T2P3M & + , T2M3P , SP0 , SP1P , SP1M , SP1P2P, SP1P2M & + , SP1P3P, SP1P3M, SP1M2P, SP1M2M, SP1M3P, SP1M3M & + , CF0 , CP0 , CF1 , CP1 , CF2 , CP2 & + , CF3 , CP3 , Q2PD0 , Q2PD1 , Q2PD2P, Q2PD3M & + , Q2MD0 , Q2MD1 , Q2MD2M, Q2MD3P ,AUX00 , AUX01 & + , AUX02 , AUX03 , AUX04 , AUX05 , SEUIL & + , AUX06 , AUX07 , AUX08 , AUX09 , AUX10 , FSEUIL + + NT = NTH + NF = NK + LBUF = 500 + SEUIL = 0. + RAISF = XFR + + DO IK = 1,NK + FREQ(IK) = FR1*RAISF**(IK-1) + ENDDO + + DO ITH = 1,NTH + DO IK = 1,NK + ! F is the E(f,theta) spectrum ... + F(ITH,IK) = DBLE(A(ITH,IK)*SIG(IK))*DBLE(TPI)/DBLE(CG(IK)) + ENDDO + ENDDO + ! CALL INSNLGQM + ! it returns: F_POIN , T_POIN , F_COEF , F_PROJ, TB_SCA , K_IF1, K_1P, k_1M , K_IF2 + ! K_IF3, K_1P2P , K_1P3M , K_1P2M , K_1P3P , K_1M2P , K_1M3M , K_1M2M + ! K_1M3P , TB_V14 , TB_FAC , TB_V24 , TB_V34 , TB_TMP , TB_TPM , IDCONF, NCONF + !======================================================================= + ! COMPUTES THE GENERALIZED MIN AND MAX FREQUENCIES : INSTEAD OF GOING + ! FROM 1 TO NF IN FREQ(JF) FOR THE MAIN FREQUENCY, IT GOES FROM JFMIN + ! TO JFMAX + ! JFMIN IS GIVEN BY Fmin=FREQ(1) /Gamma_min + ! JFMAX IS GIVEN BY Fmax=FREQ(NF)*Gamma_max + ! TESTS HAVE SHOWN THAT IT CAN BE ASSUMED Gamma_min=1. (JFMIN=1) AND + ! 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)) + ! + !======================================================================= + ! COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT + ! CALCULATED). THE THRESHOLD IS SET WITHIN 0 AND 1. + ! This was commented by FA + !======================================================================= + ! AUX00=0.0D0 + ! DO JF=1,NF + ! DO JT=1,NT + ! IF (F(JT,JF).GT.AUX00) AUX00=F(JT,JF) + ! ENDDO + ! ENDDO + ! FSEUIL=AUX00*SEUIL + + TSTOT = 0. + TSDER = 0. + !======================================================================= + ACCMAX=0. + DO JF=JFMIN,JFMAX + SUME=SUM(F(:,JF))*DTH + SATVAL(JF) = SUME*FREQ(JF)**5 + ACCVAL = SUME*FREQ(JF)**4 + IF (ACCVAL.GT.ACCMAX) ACCMAX=ACCVAL + END DO + + + ! ================================================== + ! STARTS LOOP 1 OVER THE SELECTED CONFIGURATIONS + ! ================================================== + DO ICONF=1,NCONF + ! ---------selected configuration characteristics + JF1 =IDCONF(ICONF,1) + JT1 =IDCONF(ICONF,2) + IQ_OM2=IDCONF(ICONF,3) + ! + ! ---------Recovers V1**4=(f1/f0)**4 + V1_4 =TB_V14(JF1) + ! ---------Recovers the shift of the frequency index on f1 + IXF1 =K_IF1(JF1) + ! ---------Recovers the direction indexes for Delat1 + KT1P =K_1P(JT1,JF1) + KT1M =K_1M(JT1,JF1) + ! ---------Recovers V2**4=(f2/f0)**4 and V3**4=(f3/f0)**4 + V2_4 =TB_V24(IQ_OM2,JT1,JF1) + V3_4 =TB_V34(IQ_OM2,JT1,JF1) + ! ---------Recovers the frequency indexes shift on f2 and f3 + IXF2 =K_IF2 (IQ_OM2,JT1,JF1) + IXF3 =K_IF3 (IQ_OM2,JT1,JF1) + ! ---------Recovers the direction indexes shift + KT1P2P=K_1P2P(IQ_OM2,JT1,JF1) + KT1P2M=K_1P2M(IQ_OM2,JT1,JF1) + KT1P3P=K_1P3P(IQ_OM2,JT1,JF1) + KT1P3M=K_1P3M(IQ_OM2,JT1,JF1) + KT1M2P=K_1M2P(IQ_OM2,JT1,JF1) + KT1M2M=K_1M2M(IQ_OM2,JT1,JF1) + KT1M3P=K_1M3P(IQ_OM2,JT1,JF1) + KT1M3M=K_1M3M(IQ_OM2,JT1,JF1) + ! ---------Recovers the coupling coefficients + T2P3M =TB_TPM(IQ_OM2,JT1,JF1) + T2M3P =TB_TMP(IQ_OM2,JT1,JF1) + ! ---------Recovers the multiplicative factor of QNL4 + FACTOR=TB_FAC(IQ_OM2,JT1,JF1) + + ! = = = = = = = = = = = = = = = = = = = = = = = = = + ! STARTS LOOP 2 OVER THE SPECTRUM FREQUENCIES + ! = = = = = = = = = = = = = = = = = = = = = = = = = + DO JF=JFMIN,JFMAX + IF (SATVAL(JF).GT.GQTHRSAT) THEN + ! + !.........Recovers the coefficient for the coupling factor + !.........Computes the coupling coefficients for the case +Delta1 (SIG=1) + SCAL_T=TB_SCA(LBUF+JF)*FACTOR + T_2P3M=T2P3M*SCAL_T + T_2M3P=T2M3P*SCAL_T + ! + !.........Frequency indexes and coefficients + JFM0=F_POIN(JF+LBUF) + CF0 =F_COEF(JF+LBUF) + CP0 =F_PROJ(JF+LBUF) + JFM1=F_POIN(JF+IXF1) + CF1 =F_COEF(JF+IXF1) + CP1 =F_PROJ(JF+IXF1) + JFM2=F_POIN(JF+IXF2) + CF2 =F_COEF(JF+IXF2) + CP2 =F_PROJ(JF+IXF2) + JFM3=F_POIN(JF+IXF3) + CF3 =F_COEF(JF+IXF3) + CP3 =F_PROJ(JF+IXF3) + ! + ! ------------------------------------------------- + ! STARTS LOOP 3 OVER THE SPECTRUM DIRECTIONS + ! ------------------------------------------------- + DO JT=1,NT + ! + !...........Direction indexes + ! direct config (+delta1) (sig =1) + JT1P =T_POIN(JT+KT1P) + JT1P2P=T_POIN(JT+KT1P2P) + JT1P2M=T_POIN(JT+KT1P2M) + JT1P3P=T_POIN(JT+KT1P3P) + JT1P3M=T_POIN(JT+KT1P3M) + ! image config (-delta1) + JT1M =T_POIN(JT+KT1M) + JT1M2P=T_POIN(JT+KT1M2P) + JT1M2M=T_POIN(JT+KT1M2M) + JT1M3P=T_POIN(JT+KT1M3P) + JT1M3M=T_POIN(JT+KT1M3M) + ! + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! STARTS LOOP 4 OVER THE MESH NODES + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + ! + SP0=F(JT,JFM0)*CF0 + ! + ! IF (SP0.GT.FSEUIL) THEN + ! + ! Config. +Delta1 (SIG=1) + ! ======================= + !...............Computes the spectrum values in 1, 2, 3 + SP1P =F(JT1P ,JFM1)*CF1 + SP1P2P=F(JT1P2P,JFM2)*CF2 + SP1P3M=F(JT1P3M,JFM3)*CF3 + SP1P2M=F(JT1P2M,JFM2)*CF2 + SP1P3P=F(JT1P3P,JFM3)*CF3 + ! + !...............Computes auxiliary products and variables + AUX01=SP0*V1_4+SP1P + AUX02=SP0*SP1P + AUX03=SP1P2P*SP1P3M + AUX04=SP1P2P*V3_4+SP1P3M*V2_4 + AUX05=SP1P2M*SP1P3P + AUX06=SP1P2M*V3_4+SP1P3P*V2_4 + AUX07=AUX02*V3_4 + AUX08=AUX02*V2_4 + ! + !...............Computes the components of the transfer term + S_2P3M=AUX03*AUX01-AUX02*AUX04 + S_2M3P=AUX05*AUX01-AUX02*AUX06 + Q_2P3M=T_2P3M*S_2P3M + Q_2M3P=T_2M3P*S_2M3P + AUX00 =Q_2P3M+Q_2M3P + ! + !...............Computes the components of the derived terms (dQ/dF) + Q2PD0 =T_2P3M*(AUX03*V1_4 - SP1P*AUX04)*CF0 + Q2PD1 =T_2P3M*(AUX03 - SP0 *AUX04)*CF1 + Q2PD2P=T_2P3M*(AUX01*SP1P3M - AUX07 )*CF2 + Q2PD3M=T_2P3M*(AUX01*SP1P2P - AUX08 )*CF3 + Q2MD0 =T_2M3P*(AUX05*V1_4 - SP1P*AUX06)*CF0 + Q2MD1 =T_2M3P*(AUX03 - SP0 *AUX06)*CF1 + Q2MD2M=T_2M3P*(AUX01*SP1P3P - AUX07 )*CF2 + Q2MD3P=T_2M3P*(AUX01*SP1P2M - AUX08 )*CF3 + AUX09=Q2PD0+Q2MD0 + AUX10=Q2PD1+Q2MD1 + ! + !...............Sum of Qnl4 term in the table TSTOT + TSTOT(JT,JFM0 )=TSTOT(JT,JFM0 )+AUX00 *CP0 + TSTOT(JT1P,JFM1 )=TSTOT(JT1P,JFM1 )+AUX00 *CP1 + TSTOT(JT1P2P,JFM2)=TSTOT(JT1P2P,JFM2)-Q_2P3M*CP2 + TSTOT(JT1P2M,JFM2)=TSTOT(JT1P2M,JFM2)-Q_2M3P*CP2 + TSTOT(JT1P3M,JFM3)=TSTOT(JT1P3M,JFM3)-Q_2P3M*CP3 + TSTOT(JT1P3P,JFM3)=TSTOT(JT1P3P,JFM3)-Q_2M3P*CP3 + ! + !...............Sum of the term dQnl4/dF in the table TSDER + TSDER(JT,JFM0)=TSDER(JT,JFM0)+AUX09 *CP0 + TSDER(JT1P,JFM1)=TSDER(JT1P,JFM1)+AUX10 *CP1 + TSDER(JT1P2P,JFM2)=TSDER(JT1P2P,JFM2)-Q2PD2P*CP2 + TSDER(JT1P2M,JFM2)=TSDER(JT1P2M,JFM2)-Q2MD2M*CP2 + TSDER(JT1P3M,JFM3)=TSDER(JT1P3M,JFM3)-Q2PD3M*CP3 + TSDER(JT1P3P,JFM3)=TSDER(JT1P3P,JFM3)-Q2MD3P*CP3 +#ifdef W3_TGQM + ! Test output to set up triplet method ... + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT, JFM0,AUX00 *CP0, F(JT,JFM0),TSTOT(JT ,JFM0) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1P, JFM1,AUX00 *CP1, F(JT1P,JFM1),TSTOT(JT1P,JFM1) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1P2P,JFM2,-Q_2P3M*CP2,F(JT1P2P,JFM2),TSTOT(JT1P2P,JFM2) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1P2M,JFM2,-Q_2M3P*CP2,F(JT1P2M,JFM2),TSTOT(JT1P2M,JFM2) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1P3M,JFM2,-Q_2P3M*CP3,F(JT1P3M,JFM3),TSTOT(JT1P3M,JFM3) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1P3P,JFM2,-Q_2M3P*CP3,F(JT1P3P,JFM3),TSTOT(JT1P3P,JFM3) + 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) + 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) + WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P2M,JFM2,-Q_2M3P*CP2,F(JT1P2M,JFM2),TSTOT(JT1P2M,JFM2) + WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P3M,JFM2,-Q_2P3M*CP3,F(JT1P3M,JFM3),TSTOT(JT1P3M,JFM3) + WRITE(995,'(5I3,3E12.3)') ICONF,JF,JT,JT1P3P,JFM2,-Q_2M3P*CP3,F(JT1P3P,JFM3),TSTOT(JT1P3P,JFM3) +#endif + ! + ! Config. -Delta1 (SIG=-1) + ! ======================== + !...............Computes the spectrum values in 1, 2, 3 + SP1M =F(JT1M ,JFM1)*CF1 + SP1M2P=F(JT1M2P,JFM2)*CF2 + SP1M3M=F(JT1M3M,JFM3)*CF3 + SP1M2M=F(JT1M2M,JFM2)*CF2 + SP1M3P=F(JT1M3P,JFM3)*CF3 + ! + !...............Computes auxiliary products and variables + AUX01=SP0*V1_4+SP1M + AUX02=SP0*SP1M + AUX03=SP1M2P*SP1M3M + AUX04=SP1M2P*V3_4+SP1M3M*V2_4 + AUX05=SP1M2M*SP1M3P + AUX06=SP1M2M*V3_4+SP1M3P*V2_4 + AUX07=AUX02*V3_4 + AUX08=AUX02*V2_4 + ! + !...............Computes the transfer term components + S_2P3M=AUX03*AUX01-AUX02*AUX04 + S_2M3P=AUX05*AUX01-AUX02*AUX06 + Q_2P3M=T_2M3P*S_2P3M + Q_2M3P=T_2P3M*S_2M3P + AUX00 =Q_2P3M+Q_2M3P ! Same as in +Delta1, can be commented out + ! + !...............Computes the derived terms components (dQ/dF) + Q2PD0 =T_2P3M*(AUX03*V1_4 - SP1M*AUX04)*CF0 + Q2PD1 =T_2P3M*(AUX03 - SP0 *AUX04)*CF1 + Q2PD2P=T_2P3M*(AUX01*SP1M3M - AUX07 )*CF2 + Q2PD3M=T_2P3M*(AUX01*SP1M2P - AUX08 )*CF3 + Q2MD0 =T_2M3P*(AUX05*V1_4 - SP1M*AUX06)*CF0 + Q2MD1 =T_2M3P*(AUX03 - SP0 *AUX06)*CF1 + Q2MD2M=T_2M3P*(AUX01*SP1M3P - AUX07 )*CF2 + Q2MD3P=T_2M3P*(AUX01*SP1M2M - AUX08 )*CF3 + AUX09=Q2PD0+Q2MD0 + AUX10=Q2PD1+Q2MD1 + ! + !...............Sum of Qnl4 term in the table TSTOT + TSTOT(JT ,JFM0)=TSTOT(JT ,JFM0)+AUX00 *CP0 + TSTOT(JT1M ,JFM1)=TSTOT(JT1M ,JFM1)+AUX00 *CP1 + TSTOT(JT1M2P,JFM2)=TSTOT(JT1M2P,JFM2)-Q_2P3M*CP2 + TSTOT(JT1M2M,JFM2)=TSTOT(JT1M2M,JFM2)-Q_2M3P*CP2 + TSTOT(JT1M3M,JFM3)=TSTOT(JT1M3M,JFM3)-Q_2P3M*CP3 + TSTOT(JT1M3P,JFM3)=TSTOT(JT1M3P,JFM3)-Q_2M3P*CP3 + ! + !...............Sum of the term dQnl4/dF in the table TSDER + TSDER(JT ,JFM0)=TSDER(JT ,JFM0)+AUX09 *CP0 + TSDER(JT1M ,JFM1)=TSDER(JT1M ,JFM1)+AUX10 *CP1 + TSDER(JT1M2P,JFM2)=TSDER(JT1M2P,JFM2)-Q2PD2P*CP2 + TSDER(JT1M2M,JFM2)=TSDER(JT1M2M,JFM2)-Q2MD2M*CP2 + TSDER(JT1M3M,JFM3)=TSDER(JT1M3M,JFM3)-Q2PD3M*CP3 + TSDER(JT1M3P,JFM3)=TSDER(JT1M3P,JFM3)-Q2MD3P*CP3 + ! +#ifdef W3_TGQM + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT, JFM0,AUX00 *CP0, F(JT,JFM0),TSTOT(JT ,JFM0) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1M, JFM1,AUX00 *CP1, F(JT1M,JFM1),TSTOT(JT1M,JFM1) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1M2P,JFM2,-Q_2P3M*CP2,F(JT1M2P,JFM2),TSTOT(JT1M2P,JFM2) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1M2M,JFM2,-Q_2M3P*CP2,F(JT1M2M,JFM2),TSTOT(JT1M2M,JFM2) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1M3M,JFM2,-Q_2P3M*CP3,F(JT1M3M,JFM3),TSTOT(JT1M3M,JFM3) + WRITE(994,'(5I3,3E12.3)') ICONF,JF,JT,JT1M3P,JFM2,-Q_2M3P*CP3,F(JT1M3P,JFM3),TSTOT(JT1M3P,JFM3) +#endif + ! + ! ENDIF ! this was the test on SEUIL + ! + ENDDO + ! ------------------------------------------------- + ! END OF LOOP 3 OVER THE SPECTRUM DIRECTIONS + ! ------------------------------------------------- + ! + ENDIF ! End of test on saturation level + ENDDO + ! = = = = = = = = = = = = = = = = = = = = = = = = = + ! END OF LOOP 2 OVER THE SPECTRUM FREQUENCIES + ! = = = = = = = = = = = = = = = = = = = = = = = = = + ! + ENDDO + ! ================================================== + ! END OF LOOP 1 OVER THE SELECTED CONFIGURATIONS + ! ================================================== + ! Applying WAM DEPTH SCALING ! to be added later ... + ! CALL q_dscale(F,WN,SIG,DTH,NK,NTH,DEPTH,q_dfac) + q_dfac=1 + + ! Amplification inspired by Lavrenov 2001, eq 10. + AMPFAC=GQAMP(4)*MIN(MAX(ACCMAX/GQAMP(2),1.)**GQAMP(1),GQAMP(3)) + !WRITE(991,*) ACCMAX,q_dfac,AMPFAC,GQAMP(1:3),SATVAL(10),SATVAL(30) + + ! Replacing Double Precision with Simple Real and scaling + TSTOTn = TSTOT*q_dfac*AMPFAC + TSDERn = TSDER*q_dfac*AMPFAC + + + ! Converting Snl(theta,f) to Snl(theta,k)/sigma + DO ITH = 1,NT + DO IK = 1,NF + TSTOTn(ITH,IK) = TSTOTn(ITH,IK)*CG(IK)/(TPI*SIG(IK)) + ENDDO + ENDDO + !CLOSE(994) + !STOP + END SUBROUTINE W3SNLGQM + + !/ ------------------------------------------------------------------- / + FUNCTION COUPLE(XK1 ,YK1 ,XK2 ,YK2 ,XK3 ,YK3 ,XK4 ,YK4) + !/ + !/ +-----------------------------------+ + !/ | WAVEWATCH III NOAA/NCEP | + !/ | M. Benoit & E. Gagnaire-Renou | + !/ | Last update : 20-Nov-2022 | + !/ +-----------------------------------+ + !/ + !/ 19-Nov-2022 : Transfer from TOMAWAC code ( version 7.xx ) + !/ + ! 1. Purpose : + ! + ! Computes the 4-wave coupling coefficient used in Snl4 + ! + ! 2. Method : + ! + ! Uses theoretical expression by Webb (1978) + ! + ! 3. Parameters : + ! + ! Parameter list + ! ---------------------------------------------------------------- + ! XK1 Real I x component of k1 wavenumber ... + ! ---------------------------------------------------------------- + ! + ! 5. Called by : + ! + ! Name Type Module Description + ! ---------------------------------------------------------------- + ! INNSLGQM Subr. W3SNL2 Prepares source term integration. + ! ---------------------------------------------------------------- + ! + ! 6. Error messages : + ! + ! None. + ! + ! 10. Source code : + ! + !/ ------------------------------------------------------------------- / + USE CONSTANTS, ONLY: GRAV + ! + DOUBLE PRECISION, INTENT(IN) :: XK1 , YK1 , XK2 , YK2 + DOUBLE PRECISION, INTENT(IN) :: XK3 , YK3 + DOUBLE PRECISION, INTENT(IN) :: XK4 , YK4 + DOUBLE PRECISION COUPLE + ! + !.....LOCAL VARIABLES + ! """""""""""""""""" + DOUBLE PRECISION RK1 , RK2 , RK3 , RK4 , WK1 , WK2 + DOUBLE PRECISION WK3 , WK4 , S12 , S13 , S14 , S23 + DOUBLE PRECISION S24 , S34 , W1P2 , Q12 , W1M3 , Q13 + DOUBLE PRECISION W1M4 , Q14 , DDD , COEF , DENO13, NUME13 + DOUBLE PRECISION DENO14, NUME14, ZERO, PI + + ! + PI = ACOS(-1.) + COEF=PI*GRAV*GRAV/4.D0 + ZERO=1.D-10 + ! + RK1=SQRT(XK1*XK1+YK1*YK1) + RK2=SQRT(XK2*XK2+YK2*YK2) + RK3=SQRT(XK3*XK3+YK3*YK3) + RK4=SQRT(XK4*XK4+YK4*YK4) + ! + WK1=SQRT(RK1) + WK2=SQRT(RK2) + WK3=SQRT(RK3) + WK4=SQRT(RK4) + ! + S12=XK1*XK2+YK1*YK2 + S13=XK1*XK3+YK1*YK3 + S14=XK1*XK4+YK1*YK4 + S23=XK2*XK3+YK2*YK3 + S24=XK2*XK4+YK2*YK4 + S34=XK3*XK4+YK3*YK4 + ! + W1P2=SQRT((XK1+XK2)*(XK1+XK2)+(YK1+YK2)*(YK1+YK2)) + W1M3=SQRT((XK1-XK3)*(XK1-XK3)+(YK1-YK3)*(YK1-YK3)) + W1M4=SQRT((XK1-XK4)*(XK1-XK4)+(YK1-YK4)*(YK1-YK4)) + Q12=(WK1+WK2)*(WK1+WK2) + Q13=(WK1-WK3)*(WK1-WK3) + Q14=(WK1-WK4)*(WK1-WK4) + ! + !.....COMPUTES THE D COEFFICIENT OF WEBB (1978) + ! """""""""""""""""""""""""""""""""""""" + DDD=2.00D0*Q12*(RK1*RK2-S12)*(RK3*RK4-S34)/(W1P2-Q12) & + +0.50D0*(S12*S34+S13*S24+S14*S23) & + +0.25D0*(S13+S24)*Q13*Q13 & + -0.25D0*(S12+S34)*Q12*Q12 & + +0.25D0*(S14+S23)*Q14*Q14 & + +2.50D0*RK1*RK2*RK3*RK4 & + +Q12*Q13*Q14*(RK1+RK2+RK3+RK4) + + DENO13=W1M3-Q13 + NUME13=2.00D0*Q13*(RK1*RK3+S13)*(RK2*RK4+S24) + IF (ABS(DENO13).LT.ZERO) THEN + IF (ABS(NUME13).LT.ZERO) THEN + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3) 0/0 !' + ELSE + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3) inifinte value' + ENDIF + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-3) term not used' + ELSE + DDD=DDD+NUME13/DENO13 + ENDIF + DENO14=W1M4-Q14 + NUME14=2.00D0*Q14*(RK1*RK4+S14)*(RK2*RK3+S23) + IF (ABS(DENO14).LT.ZERO) THEN + IF (ABS(NUME14).LT.ZERO) THEN + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4) 0/0 !' + ELSE + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4) inifinte value' + ENDIF + WRITE(*,*) 'W3SNL2 error for coupling coefficient : (1-4) term not used' + ELSE + DDD=DDD+NUME14/DENO14 + ENDIF + + COUPLE=COEF*DDD*DDD/(WK1*WK2*WK3*WK4) + ! RETURN + END FUNCTION COUPLE + + !/ ------------------------------------------------------------------- / + SUBROUTINE GAULEG (W_LEG ,X_LEG ,NPOIN) + !/ ------------------------------------------------------------------- / + !.....VARIABLES IN ARGUMENT + ! """""""""""""""""""" + INTEGER , INTENT(IN) :: NPOIN + DOUBLE PRECISION ,INTENT(INOUT) :: W_LEG(NPOIN) , X_LEG(NPOIN) + ! + !.....LOCAL VARIABLES + ! """"""""""""""""" + INTEGER I, M, J + DOUBLE PRECISION EPS, Z, P1, P2, P3, PP, Z1, PI + PARAMETER (EPS=3.D-14) + ! + PI = ACOS(-1.) + M=(NPOIN+1)/2 + DO I=1,M + Z=COS(PI*(DBLE(I)-0.25D0)/(DBLE(NPOIN)+0.5D0)) +1 CONTINUE + P1=1.0D0 + P2=0.0D0 + DO J=1,NPOIN + P3=P2 + P2=P1 + P1=((2.D0*DBLE(J)-1.D0)*Z*P2-(DBLE(J)-1.D0)*P3)/DBLE(J) + ENDDO + PP=DBLE(NPOIN)*(Z*P1-P2)/(Z*Z-1.D0) + Z1=Z + Z=Z-P1/PP + IF (ABS(Z-Z1).GT.EPS) GOTO 1 + X_LEG(I)=-Z + X_LEG(NPOIN+1-I)=Z + W_LEG(I)=2.D0/((1.D0-Z**2)*PP**2) + W_LEG(NPOIN+1-I)=W_LEG(I) + ENDDO + END SUBROUTINE GAULEG + + !/ ------------------------------------------------------------------- / + SUBROUTINE F1F1F1(F1SF,NF1,IQ_OM1) + ! TOMAWAC V6P3 15/06/2011 + !*********************************************************************** + ! + !brief SUBROUTINE CALLED BY PRENL3 + !+ COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1 + !+ INDICATOR + ! + !history E. GAGNAIRE-RENOU + !+ 04/2011 + !+ V6P1 + !+ CREATED + ! + !history G.MATTAROLO (EDF - LNHE) + !+ 15/06/2011 + !+ V6P1 + !+ Translation of French names of the variables in argument + ! + !history E. GAGNAIRE-RENOU + !+ 12/03/2013 + !+ V6P3 + !+ Better formatted: WRITE(LU,*), etc. + !/ ------------------------------------------------------------------- / + IMPLICIT NONE + INTEGER, INTENT(IN) :: IQ_OM1 + INTEGER, INTENT(INOUT) :: NF1 + DOUBLE PRECISION, INTENT(INOUT) :: F1SF(*) + ! + INTEGER I,M + DOUBLE PRECISION RAISON + ! + IF(IQ_OM1.EQ.1) THEN + IF(NF1.NE.14) THEN + WRITE(*,*) '#1 Incorrect value for NF1',NF1 + ENDIF + F1SF( 1)=0.30D0 + F1SF( 2)=0.40D0 + F1SF( 3)=0.50D0 + F1SF( 4)=0.60D0 + F1SF( 5)=0.70D0 + F1SF( 6)=0.80D0 + F1SF( 7)=0.90D0 + F1SF( 8)=1.00D0 + F1SF( 9)=1.11D0 + F1SF(10)=1.25D0 + F1SF(11)=1.42D0 + F1SF(12)=1.67D0 + F1SF(13)=2.00D0 + F1SF(14)=2.50D0 + F1SF(15)=3.30D0 + ELSEIF(IQ_OM1.EQ.2) THEN + IF (NF1.NE.26) THEN + WRITE(*,*) '#2 Incorrect value for NF1', NF1 + ENDIF + F1SF( 1)=0.32D0 + F1SF( 2)=0.35D0 + F1SF( 3)=0.39D0 + F1SF( 4)=0.44D0 + F1SF( 5)=0.50D0 + F1SF( 6)=0.56D0 + F1SF( 7)=0.63D0 + F1SF( 8)=0.70D0 + F1SF( 9)=0.78D0 + F1SF(10)=0.86D0 + F1SF(11)=0.92D0 + F1SF(12)=0.97D0 + F1SF(13)=1.00D0 + F1SF(14)=1.03D0 + F1SF(15)=1.08D0 + F1SF(16)=1.13D0 + F1SF(17)=1.20D0 + F1SF(18)=1.28D0 + F1SF(19)=1.37D0 + F1SF(20)=1.48D0 + F1SF(21)=1.50D0 + F1SF(22)=1.65D0 + F1SF(23)=1.85D0 + F1SF(24)=2.10D0 + F1SF(25)=2.40D0 + F1SF(26)=2.70D0 + F1SF(27)=3.20D0 + ELSEIF(IQ_OM1.EQ.3) THEN + IF(NF1.NE.11) THEN + WRITE(*,*) 'Incorrect value for NF1', NF1 + ENDIF + F1SF( 1)=0.30D0 + F1SF( 2)=0.48D0 + F1SF( 3)=0.64D0 + F1SF( 4)=0.78D0 + F1SF( 5)=0.90D0 + F1SF( 6)=1.00D0 + F1SF( 7)=1.12D0 + F1SF( 8)=1.28D0 + F1SF( 9)=1.50D0 + F1SF(10)=1.80D0 + F1SF(11)=2.40D0 + F1SF(12)=3.40D0 + ELSEIF(IQ_OM1.EQ.4) THEN + IF(NF1.NE.40) THEN + WRITE(*,*) 'Incorrect value for NF1', NF1 + ENDIF + NF1=20 + M=10 + RAISON=9.D0**(1.D0/DBLE(NF1)) + F1SF(M+1)=1.0D0/3.0D0 + NF1=2*M+NF1 + DO I=M+2,NF1+1 + F1SF(I)=F1SF(I-1)*RAISON + ENDDO + DO I=M,1,-1 + F1SF(I)=F1SF(I+1)/RAISON + ENDDO + ELSEIF(IQ_OM1.EQ.5) THEN + RAISON=9.D0**(1.D0/DBLE(NF1)) + F1SF(1)=1.D0/3.D0 + DO I=2,NF1+1 + F1SF(I)=F1SF(I-1)*RAISON + ENDDO + ELSEIF(IQ_OM1.EQ.6) THEN + RAISON=(3.D0-1.D0/3.D0)/DBLE(NF1) + F1SF(1)=1.D0/3.D0 + DO I=2,NF1+1 + F1SF(I)=F1SF(I-1)+RAISON + ENDDO + ELSEIF(IQ_OM1.EQ.7) THEN + IF(NF1.NE.20) THEN + WRITE(*,*) 'Incorrect value for NF1', NF1 + ENDIF + F1SF( 1)=1.D0/3.D0 + F1SF( 2)=0.40D0 + F1SF( 3)=0.46D0 + F1SF( 4)=0.52D0 + F1SF( 5)=0.60D0 + F1SF( 6)=0.70D0 + F1SF( 7)=0.79D0 + F1SF( 8)=0.86D0 + F1SF( 9)=0.92D0 + F1SF(10)=0.97D0 + F1SF(11)=1.00D0 + F1SF(12)=1.04D0 + F1SF(13)=1.10D0 + F1SF(14)=1.18D0 + F1SF(15)=1.28D0 + F1SF(16)=1.42D0 + F1SF(17)=1.60D0 + F1SF(18)=1.84D0 + F1SF(19)=2.14D0 + F1SF(20)=2.52D0 + F1SF(21)=3.00D0 + ENDIF + ! + END SUBROUTINE F1F1F1 + !/ ------------------------------------------------------------------- / + SUBROUTINE INSNLGQM + !/ + !/ +-----------------------------------+ + !/ | WAVEWATCH III NOAA/NCEP | + !/ | E. Gagnaire-Renou & | + !/ | M. Benoit | + !/ | S. Mostafa Siadatamousavi | + !/ | M. Beyramzadeh | + !/ | FORTRAN 90 | + !/ | Last update : 20-Nov-2022 | + !/ +-----------------------------------+ + !/ + !/ 20-Nov-2022 : Merging with NL2 in WW3. ( version 7.00 ) + !/ + ! 1. Purpose : + ! + ! Preprocessing for nonlinear interactions (Xnl). + ! + ! 2. Method : + ! + ! See Xnl documentation. + ! + ! 3. Parameters : + ! + ! 4. Subroutines used : + ! + ! Name Type Module Description + ! ---------------------------------------------------------------- + ! STRACE Subr. W3SERVMD Subroutine tracing. + ! Subr. GAULEG Gauss-Legendre weights + ! xnl_init Subr. m_constants Xnl initialization routine. + ! ---------------------------------------------------------------- + ! + ! 5. Called by : + ! + ! Name Type Module Description + ! ---------------------------------------------------------------- + ! W3IOGR Subr. W3IOGRMD Model definition file management. + ! ---------------------------------------------------------------- + ! + ! 6. Error messages : + ! + ! 7. Remarks : + ! + ! 8. Structure : + ! + ! - See source code. + ! + ! 9. Switches : + ! + ! !/S Enable subroutine tracing. + ! + ! 10. Source code : + ! + !/ ------------------------------------------------------------------- / + USE CONSTANTS, ONLY: GRAV + USE W3GDATMD, ONLY: NK , NTH , XFR , FR1, GQNF1, GQNT1, GQNQ_OM2, NLTAIL, GQTHRCOU + +#ifdef W3_S + CALL STRACE (IENT, 'INSNLGQM') +#endif + !.....LOCAL VARIABLES + INTEGER JF , JT , JF1 , JT1 , NF1P1 , IAUX , NT , NF , IK + INTEGER IQ_TE1 , IQ_OM2 , LBUF , DIMBUF , IQ_OM1 , NQ_TE1 , NCONFM + + DOUBLE PRECISION EPSI_A, AUX , CCC , DENO , AAA , DP2SG , TAILF + DOUBLE PRECISION V1 , V1_4 , DV1 , DTETAR , ELIM , RAISF + DOUBLE PRECISION V2 , V2_4 , V3 , V3_4 + DOUBLE PRECISION W2 , W2_M , W2_1 , W_MIL , W_RAD + DOUBLE PRECISION RK0 , XK0 , YK0 , RK1 , XK1 , YK1 + DOUBLE PRECISION RK2 , XK2P , YK2P , XK2M , YK2M + DOUBLE PRECISION RK3 , XK3P , YK3P , XK3M , YK3M + DOUBLE PRECISION D01P , C_D01P, S_D01P, D0AP , C_D0AP, S_D0AP + DOUBLE PRECISION GA2P , C_GA2P, S_GA2P, GA3P , C_GA3P, S_GA3P, TWOPI, PI, SEUIL1 , SEUIL2 , SEUIL + ! + !.....Variables related to the Gaussian quadratures + DOUBLE PRECISION W_CHE_TE1, W_CHE_OM2, C_LEG_OM2 + ! + !.....Variables related to the configuration selection + DOUBLE PRECISION TEST1 , TEST2 + DOUBLE PRECISION :: FREQ(NK) + DOUBLE PRECISION, ALLOCATABLE :: F1SF(:) , X_CHE_TE1(:) , X_CHE_OM2(:) , X_LEG_OM2(:) , W_LEG_OM2(:) & + , MAXCLA(:) + + PI = Acos(-1.) + LBUF = 500 + DIMBUF = 2*LBUF+200 + TWOPI = 2.*PI + ! + ! Defines some threshold values for filtering (See Gagnaire-Renou Thesis, p 52) + ! + SEUIL1 = 1E10 + SEUIL2 = GQTHRCOU + + IF(GQNF1.EQ.14) IQ_OM1=1 + IF(GQNF1.EQ.26) IQ_OM1=2 + IF(GQNF1.EQ.11) IQ_OM1=3 + IF(GQNF1.EQ.40) IQ_OM1=4 + IF(GQNF1.EQ.11) IQ_OM1=3 + IF(GQNF1.EQ.40) IQ_OM1=4 + IF(GQNF1.EQ.20) IQ_OM1=7 + ! + ! Note by FA: not sure what the 5 and 6 cases correspond to + ! + NQ_TE1 = GQNT1/2 + NCONFM = GQNF1*GQNT1*GQNQ_OM2 + + RAISF = XFR + NT = NTH + NF = NK + DTETAR = TWOPI/DBLE(NT) + + DO IK = 1,NK + FREQ(IK) = FR1*RAISF**(IK-1) + ENDDO + + TAILF = -NLTAIL + + !===============ALLOCATE MATRICES============================================= + if (Allocated(K_IF2) ) then + deallocate(K_IF2) + endif + ALLOCATE(K_IF2(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_IF3) ) then + deallocate(K_IF3) + endif + ALLOCATE(K_IF3(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1P2P) ) then + deallocate(K_1P2P) + endif + ALLOCATE(K_1P2P(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1P3M) ) then + deallocate(K_1P3M) + endif + ALLOCATE(K_1P3M(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1P2M) ) then + deallocate(K_1P2M) + endif + ALLOCATE(K_1P2M(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1P3P) ) then + deallocate(K_1P3P) + endif + ALLOCATE(K_1P3P(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1M2P) ) then + deallocate(K_1M2P) + endif + ALLOCATE(K_1M2P(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1M3M) ) then + deallocate(K_1M3M) + endif + ALLOCATE(K_1M3M(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1M2M) ) then + deallocate(K_1M2M) + endif + ALLOCATE(K_1M2M(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_1M3P) ) then + deallocate(K_1M3P) + endif + ALLOCATE(K_1M3P(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(TB_V24) ) then + deallocate(TB_V24) + endif + ALLOCATE(TB_V24(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(TB_V34) ) then + deallocate(TB_V34) + endif + ALLOCATE(TB_V34(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(TB_TPM) ) then + deallocate(TB_TPM) + endif + ALLOCATE(TB_TPM(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(TB_TMP) ) then + deallocate(TB_TMP) + endif + ALLOCATE(TB_TMP(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(TB_FAC) ) then + deallocate(TB_FAC) + endif + ALLOCATE(TB_FAC(GQNQ_OM2,GQNT1,GQNF1)) + + if (Allocated(K_IF1) ) then + deallocate(K_IF1) + endif + ALLOCATE(K_IF1(GQNF1)) + + if (Allocated(K_1P) ) then + deallocate(K_1P) + endif + ALLOCATE(K_1P(GQNT1,GQNF1)) + + if (Allocated(K_1M) ) then + deallocate(K_1M) + endif + ALLOCATE(K_1M(GQNT1,GQNF1)) + + if (Allocated(TB_V14) ) then + deallocate(TB_V14) + endif + ALLOCATE(TB_V14(GQNF1)) + + if (Allocated(IDCONF) ) then + deallocate(IDCONF) + endif + ALLOCATE(IDCONF(NCONFM,3)) + + !======================================================================= + ! INITIALISATION OF AUXILIAIRY TABLES FOR SPECTRUM INTERPOLATION + !======================================================================= + if (Allocated(F_POIN) ) then + deallocate(F_POIN) + endif + ALLOCATE(F_POIN(DIMBUF)) + + if (Allocated(T_POIN) ) then + deallocate(T_POIN) + endif + ALLOCATE(T_POIN(DIMBUF)) + + if (Allocated(F_COEF) ) then + deallocate(F_COEF) + endif + ALLOCATE(F_COEF(DIMBUF)) + + if (Allocated(F_PROJ) ) then + deallocate(F_PROJ) + endif + ALLOCATE(F_PROJ(DIMBUF)) + + if (Allocated(TB_SCA) ) then + deallocate(TB_SCA) + endif + ALLOCATE(TB_SCA(DIMBUF)) + + + F_POIN(:)=0 + T_POIN(:)=0 + F_COEF(:)=0.D0 + F_PROJ(:)=0.D0 + TB_SCA(:)=0.0D0 + + DO JF=1,LBUF + F_POIN(JF)=1 + F_COEF(JF)=0.0D0 + F_PROJ(JF)=0.0D0 + ENDDO + DO JF=1,NF + IAUX=LBUF+JF + F_POIN(IAUX)=JF + F_COEF(IAUX)=1.0D0 + F_PROJ(IAUX)=1.0D0 + ENDDO + AUX=1.D0/RAISF**TAILF + DO JF=1,LBUF + IAUX=LBUF+NF+JF + F_POIN(IAUX)=NF + F_COEF(IAUX)=AUX**JF + F_PROJ(IAUX)=0.0D0 + ENDDO + ! + DO JT=LBUF,1,-1 + T_POIN(JT)=NT-MOD(LBUF-JT,NT) + ENDDO + DO JT=1,NT + T_POIN(LBUF+JT)=JT + ENDDO + DO JT=1,LBUF + T_POIN(LBUF+NT+JT)=MOD(JT-1,NT)+1 + ENDDO + !====================================================================== + ! + !======================================================================= + ! COMPUTES SCALE COEFFICIENTS FOR THE COUPLING COEFFICIENT + ! Would be easier to pass these on from W3SRCE ??? + !======================================================================= + DP2SG=TWOPI*TWOPI/GRAV + DO JF=1,LBUF + AUX=FREQ(1)/RAISF**(LBUF-JF+1) + TB_SCA(JF)=(DP2SG*AUX**2)**6/(TWOPI**3*AUX) + ENDDO + DO JF=1,NF + TB_SCA(LBUF+JF)=(DP2SG*FREQ(JF)**2)**6/(TWOPI**3*FREQ(JF)) + ENDDO + DO JF=1,LBUF + IAUX=LBUF+NF+JF + AUX=FREQ(NF)*RAISF**JF + TB_SCA(IAUX)=(DP2SG*AUX**2)**6/(TWOPI**3*AUX) + ENDDO + !======================================================================= + ! + !======================================================================= + ! COMPUTES VALUES FOR GAUSSIAN QUADRATURES + !======================================================================= + if (Allocated(X_CHE_TE1) ) then + deallocate(X_CHE_TE1) + endif + ALLOCATE(X_CHE_TE1(1:NQ_TE1),X_CHE_OM2(1:GQNQ_OM2)) + + if (Allocated(X_LEG_OM2) ) then + deallocate(X_LEG_OM2) + endif + ALLOCATE(X_LEG_OM2(1:GQNQ_OM2),W_LEG_OM2(1:GQNQ_OM2)) + ! + !.....Abscissa and weight (constant) for Gauss-Chebyshev + DO IQ_TE1=1,NQ_TE1 + X_CHE_TE1(IQ_TE1)=COS(PI*(DBLE(IQ_TE1)-0.5D0)/DBLE(NQ_TE1)) + ENDDO + W_CHE_TE1=PI/DBLE(NQ_TE1) + DO IQ_OM2=1,GQNQ_OM2 + X_CHE_OM2(IQ_OM2)=COS(PI*(DBLE(IQ_OM2)-0.5D0)/DBLE(GQNQ_OM2)) + ENDDO + W_CHE_OM2=PI/DBLE(GQNQ_OM2) + ! + !.....Abscissa et weight for Gauss-Legendre + CALL GAULEG( W_LEG_OM2 , X_LEG_OM2 , GQNQ_OM2 ) + DO IQ_OM2=1,GQNQ_OM2 + X_LEG_OM2(IQ_OM2)=0.25D0*(1.D0+X_LEG_OM2(IQ_OM2))**2 + ENDDO + !======================================================================= + ! + ! + !======================================================================= + ! COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1 INDICATOR + !======================================================================= + NF1P1=GQNF1+1 + if (Allocated(F1SF) ) then + deallocate(F1SF) + endif + ALLOCATE(F1SF(1:NF1P1)) + + CALL F1F1F1 ( F1SF , GQNF1 , IQ_OM1) + !======================================================================= + ! + ! ================================================== + ! STARTS LOOP 1 OVER THE RATIOS F1/F0 + ! ================================================== + DO JF1=1,GQNF1 + ! ---------Computes and stores v1=f1/f0 and v1**4 + V1=(F1SF(JF1+1)+F1SF(JF1))/2.D0 + K_IF1(JF1)=NINT(DBLE(LBUF)+LOG(V1)/LOG(RAISF)) + V1_4=V1**4 + TB_V14(JF1)=V1_4 + ! ---------Computes and stores dv1=df1/f0 + DV1=F1SF(JF1+1)-F1SF(JF1) + ! ---------Computes the A parameter + AAA=((1.D0+V1)**4-4.D0*(1.D0+V1_4))/(8.D0*V1**2) + ! + ! ================================================= + ! STARTS LOOP 2 OVER THE DELTA_1+ VALUES + ! ================================================= + DO JT1=1,GQNT1 + ! + !......Computes the Delta1+ values (=Theta_1-Theta_0) between 0 and Pi. + IF (JT1.LE.NQ_TE1) THEN + ! ---------First interval : X from -1 to A + IQ_TE1=JT1 + C_D01P=(-1.D0+AAA)/2.D0+(1.D0+AAA)/2.D0*X_CHE_TE1(IQ_TE1) + CCC=DV1*SQRT((AAA-C_D01P)/(1.D0-C_D01P))*W_CHE_TE1 + ELSE + ! ---------Second interval : X from A to 1 + IQ_TE1=JT1-NQ_TE1 + C_D01P=( 1.D0+AAA)/2.D0+(1.D0-AAA)/2.D0*X_CHE_TE1(IQ_TE1) + CCC=DV1*SQRT((C_D01P-AAA)/(1.D0+C_D01P))*W_CHE_TE1 + ENDIF + S_D01P=SQRT(1.D0-C_D01P*C_D01P) + D01P =ACOS(C_D01P) + K_1P(JT1,JF1)=LBUF+NINT(D01P/DTETAR) + K_1M(JT1,JF1)=LBUF-NINT(D01P/DTETAR) + ! + ! ---------Computes Epsilon_a + EPSI_A=2.D0*SQRT(1.D0+V1_4+2.D0*V1*V1*C_D01P)/(1.D0+V1)**2 + ! ---------Computes Delta_A+ and its cosinus + C_D0AP=(1.D0-V1_4+0.25D0*EPSI_A**2*(1.D0+V1)**4) & + /(EPSI_A*(1.D0+V1)**2) + S_D0AP=SQRT(1.0D0-C_D0AP*C_D0AP) + D0AP = ACOS(C_D0AP) + ! + !.......Integration over OMEGA2 depending on EPS_A + IF (EPSI_A.LT.1.D0) THEN + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !........Case of a single singularity (in OMEGA2-) + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - + W2_M=0.5D0*(1.D0-EPSI_A/2.D0) + W2_1=0.5D0 + ! + W_RAD=W2_1-W2_M + C_LEG_OM2=SQRT(W_RAD) + ! + ! ---------------------------------------------------- + !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1) + !........Case of a single singularity (in OMEGA2-) + !........Integration over OMEGA2 via GAUSS-LEGENDRE quadrature + ! ---------------------------------------------------- + DO IQ_OM2=1,GQNQ_OM2 + ! ---------Computes W2, V2, and V3 + W2=W2_M+W_RAD*X_LEG_OM2(IQ_OM2) + V2=W2*(1.D0+V1) + V2_4=V2**4 + TB_V24(IQ_OM2,JT1,JF1)=V2_4 + K_IF2 (IQ_OM2,JT1,JF1) = NINT(DBLE(LBUF) & + + LOG(V2)/LOG(RAISF)) + V3=1.D0+V1-V2 + V3_4=V3**4 + TB_V34(IQ_OM2,JT1,JF1)=V3_4 + K_IF3 (IQ_OM2,JT1,JF1) = NINT(DBLE(LBUF) & + + LOG(V3)/LOG(RAISF)) + ! ---------Computes Gamma_2+ et Gamma_3+ angles + C_GA2P=(EPSI_A**2/4.D0+W2**4-(1.D0-W2)**4)/(EPSI_A*W2*W2) + C_GA2P=MAX(MIN(C_GA2P,1.D0),-1.D0) + S_GA2P=SQRT(1.D0-C_GA2P*C_GA2P) + GA2P =ACOS(C_GA2P) + C_GA3P=(EPSI_A**2/4.D0-W2**4+(1.D0-W2)**4)/EPSI_A & + /(1.D0-W2)**2 + C_GA3P=MAX(MIN(C_GA3P,1.D0),-1.D0) + S_GA3P=SQRT(1.D0-C_GA3P*C_GA3P) + GA3P =ACOS(C_GA3P) + ! Shifting of the direction indexes - Config. +Delta1 (SIG=1) + K_1P2P(IQ_OM2,JT1,JF1)=NINT(( D0AP+GA2P)/DTETAR & + +DBLE(LBUF)) + K_1P3M(IQ_OM2,JT1,JF1)=NINT(( D0AP-GA3P)/DTETAR & + +DBLE(LBUF)) + K_1P2M(IQ_OM2,JT1,JF1)=NINT(( D0AP-GA2P)/DTETAR & + +DBLE(LBUF)) + K_1P3P(IQ_OM2,JT1,JF1)=NINT(( D0AP+GA3P)/DTETAR & + +DBLE(LBUF)) + ! Shifting of the direction indexes - Config. -Delta1 (SIG=-1) + K_1M2P(IQ_OM2,JT1,JF1)=NINT((-D0AP+GA2P)/DTETAR & + +DBLE(LBUF)) + K_1M3M(IQ_OM2,JT1,JF1)=NINT((-D0AP-GA3P)/DTETAR & + +DBLE(LBUF)) + K_1M2M(IQ_OM2,JT1,JF1)=NINT((-D0AP-GA2P)/DTETAR & + +DBLE(LBUF)) + K_1M3P(IQ_OM2,JT1,JF1)=NINT((-D0AP+GA3P)/DTETAR & + +DBLE(LBUF)) + ! + !.........Computes the coupling coefficients (only for Delta_1+ ) + RK0=1.D0 + RK1=V1*V1 + RK2=V2*V2 + RK3=(1.D0+V1-V2)**2 + XK0 = RK0 + YK0 = 0.0D0 + XK1 = RK1*C_D01P + YK1 = RK1*S_D01P + XK2P = RK2*(C_D0AP*C_GA2P-S_D0AP*S_GA2P) + YK2P = RK2*(S_D0AP*C_GA2P+C_D0AP*S_GA2P) + XK2M = RK2*(C_D0AP*C_GA2P+S_D0AP*S_GA2P) + YK2M = RK2*(S_D0AP*C_GA2P-C_D0AP*S_GA2P) + XK3P = RK3*(C_D0AP*C_GA3P-S_D0AP*S_GA3P) + YK3P = RK3*(S_D0AP*C_GA3P+C_D0AP*S_GA3P) + XK3M = RK3*(C_D0AP*C_GA3P+S_D0AP*S_GA3P) + YK3M = RK3*(S_D0AP*C_GA3P-C_D0AP*S_GA3P) + TB_TPM(IQ_OM2,JT1,JF1)=COUPLE( XK0 , YK0 , XK1 , YK1 , XK2P , YK2P , XK3M , YK3M) + TB_TMP(IQ_OM2,JT1,JF1)=COUPLE( XK0 , YK0 , XK1 , YK1 , XK2M , YK2M , XK3P , YK3P) + ! + !.........Computes the multiplicative coefficient for QNL4 + DENO=2.D0*SQRT( (0.5D0*(1.D0+EPSI_A/2.D0)-W2) & + *((W2-0.5D0)**2+0.25D0*(1.D0+EPSI_A)) & + *((W2-0.5D0)**2+0.25D0*(1.D0-EPSI_A)) ) + TB_FAC(IQ_OM2,JT1,JF1)=1.D0/(DENO*V1*W2*(1.D0-W2)) & + /(1.D0+V1)**5 * W_LEG_OM2(IQ_OM2)*C_LEG_OM2* CCC + ENDDO + ! ----------------------------------------------- + !........END OF THE LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1) + ! ----------------------------------------------- + ! + ELSE + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - + !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1) + !........Case of two singularities (in OMEGA2- and OMEGA2_1) + !........Integration over OMEGA2 via GAUSS-CHEBYSCHEV quadrature + ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - + W2_M=0.5D0*(1.D0-EPSI_A/2.D0) + W2_1=0.5D0*(1.D0-SQRT(EPSI_A-1.D0)) + ! + W_MIL=(W2_M+W2_1)/2.D0 + W_RAD=(W2_1-W2_M)/2.D0 + ! + DO IQ_OM2=1,GQNQ_OM2 + ! ---------Computes W2, V2, and V3 + W2=W_MIL+W_RAD*X_CHE_OM2(IQ_OM2) + V2=W2*(1.D0+V1) + V2_4=V2**4 + TB_V24(IQ_OM2,JT1,JF1)=V2_4 + K_IF2 (IQ_OM2,JT1,JF1)=NINT(DBLE(LBUF) & + +LOG(V2)/LOG(RAISF)) + V3=1.D0+V1-V2 + V3_4=V3**4 + TB_V34(IQ_OM2,JT1,JF1)=V3_4 + K_IF3 (IQ_OM2,JT1,JF1)=NINT(DBLE(LBUF) & + +LOG(V3)/LOG(RAISF)) + ! ---------Computes Gamma_2+ et Gamma_3+ angles + C_GA2P=(EPSI_A**2/4.D0+W2**4-(1.D0-W2)**4)/(EPSI_A*W2*W2) + C_GA2P=MAX(MIN(C_GA2P,1.D0),-1.D0) + S_GA2P=SQRT(1.D0-C_GA2P*C_GA2P) + GA2P =ACOS(C_GA2P) + C_GA3P=(EPSI_A**2/4.D0-W2**4+(1.D0-W2)**4)/EPSI_A & + /(1.D0-W2)**2 + C_GA3P=MAX(MIN(C_GA3P,1.D0),-1.D0) + S_GA3P=SQRT(1.D0-C_GA3P*C_GA3P) + GA3P =ACOS(C_GA3P) + ! Shifts the direction indexes - Config. +Delta1 (SIG=1) + K_1P2P(IQ_OM2,JT1,JF1)=NINT(( D0AP+GA2P)/DTETAR & + +DBLE(LBUF)) + K_1P3M(IQ_OM2,JT1,JF1)=NINT(( D0AP-GA3P)/DTETAR & + +DBLE(LBUF)) + K_1P2M(IQ_OM2,JT1,JF1)=NINT(( D0AP-GA2P)/DTETAR & + +DBLE(LBUF)) + K_1P3P(IQ_OM2,JT1,JF1)=NINT(( D0AP+GA3P)/DTETAR & + +DBLE(LBUF)) + ! Shifts the direction indexes - Config. -Delta1 (SIG=-1) + K_1M2P(IQ_OM2,JT1,JF1)=NINT((-D0AP+GA2P)/DTETAR & + +DBLE(LBUF)) + K_1M3M(IQ_OM2,JT1,JF1)=NINT((-D0AP-GA3P)/DTETAR & + +DBLE(LBUF)) + K_1M2M(IQ_OM2,JT1,JF1)=NINT((-D0AP-GA2P)/DTETAR & + +DBLE(LBUF)) + K_1M3P(IQ_OM2,JT1,JF1)=NINT((-D0AP+GA3P)/DTETAR & + +DBLE(LBUF)) + ! + !.........Computes the coupling coefficients (only for Delta_1+ ) + RK0=1.D0 + RK1=V1*V1 + RK2=V2*V2 + RK3=(1.D0+V1-V2)**2 + XK0 = RK0 + YK0 = 0.0D0 + XK1 = RK1*C_D01P + YK1 = RK1*S_D01P + XK2P = RK2*(C_D0AP*C_GA2P-S_D0AP*S_GA2P) + YK2P = RK2*(S_D0AP*C_GA2P+C_D0AP*S_GA2P) + XK2M = RK2*(C_D0AP*C_GA2P+S_D0AP*S_GA2P) + YK2M = RK2*(S_D0AP*C_GA2P-C_D0AP*S_GA2P) + XK3P = RK3*(C_D0AP*C_GA3P-S_D0AP*S_GA3P) + YK3P = RK3*(S_D0AP*C_GA3P+C_D0AP*S_GA3P) + XK3M = RK3*(C_D0AP*C_GA3P+S_D0AP*S_GA3P) + YK3M = RK3*(S_D0AP*C_GA3P-C_D0AP*S_GA3P) + TB_TPM(IQ_OM2,JT1,JF1)=COUPLE( XK0 , YK0 , XK1 , YK1 , XK2P , YK2P , XK3M , YK3M) + TB_TMP(IQ_OM2,JT1,JF1)=COUPLE( XK0 , YK0 , XK1 , YK1 , XK2M , YK2M , XK3P , YK3P) + ! + !.........Computes the multiplicative coefficient for QNL4 + DENO=2.D0*SQRT( (0.5D0*(1.D0+EPSI_A/2.D0)-W2) & + *((W2-0.5D0)**2+0.25D0*(1.D0+EPSI_A)) & + *(0.5D0*(1.D0+SQRT(EPSI_A-1.D0))-W2) ) + TB_FAC(IQ_OM2,JT1,JF1)=1.D0/(DENO*V1*W2*(1.D0-W2)) & + /(1.D0+V1)**5 * W_CHE_OM2* CCC + ! + ENDDO + ! ----------------------------------------------- + !........END OF LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1) + ! ----------------------------------------------- + ! + ENDIF + ENDDO + ! ================================================= + ! END OF LOOP 2 OVER THE DELTA_1+ VALUES + ! ================================================= + ! + ENDDO + ! ================================================== + ! END OF LOOP 1 OVER THE F1/F0 RATIOS + ! ================================================== + DEALLOCATE(F1SF) + DEALLOCATE(X_CHE_TE1) + DEALLOCATE(X_CHE_OM2) + DEALLOCATE(X_LEG_OM2) + DEALLOCATE(W_LEG_OM2) + + ! =========================================================== + ! POST-PROCESSING TO ELIMINATE PART OF THE CONFIGURATIONS + ! =========================================================== + ! + !.....It looks, for every value of the ratio V1, for the maximum value + !.....of FACTOR*COUPLING : it is stored in the local table NAXCLA(.) + ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + ALLOCATE(MAXCLA(1:GQNF1)) + DO JF1=1,GQNF1 + 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 + ENDDO + ENDDO + MAXCLA(JF1)=AUX + ENDDO + ! + !.....It looks for the max V1 value + ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + AUX=0.0D0 + DO JF1=1,GQNF1 + IF (MAXCLA(JF1).GT.AUX) AUX=MAXCLA(JF1) + ENDDO + TEST1=SEUIL1*AUX + ! + !.....Set to zero the coupling coefficients not used + ! """"""""""""""""""""""""""""""""""""""""""""""""""""" + NCONF=0 + DO JF1=1,GQNF1 + TEST2 =SEUIL2*MAXCLA(JF1) + DO JT1=1,GQNT1 + DO IQ_OM2=1,GQNQ_OM2 + AAA=TB_FAC(IQ_OM2,JT1,JF1)*TB_TPM(IQ_OM2,JT1,JF1) + CCC=TB_FAC(IQ_OM2,JT1,JF1)*TB_TMP(IQ_OM2,JT1,JF1) + IF ((AAA.GT.TEST1.OR.AAA.GT.TEST2).OR. & + (CCC.GT.TEST1.OR.CCC.GT.TEST2)) THEN + NCONF=NCONF+1 + IDCONF(NCONF,1)=JF1 + IDCONF(NCONF,2)=JT1 + IDCONF(NCONF,3)=IQ_OM2 + ENDIF +#ifdef W3_TGQM + WRITE(993,*) NCONF,JF1,JT1,IQ_OM2,AAA,CCC,(AAA.GT.TEST1.OR.AAA.GT.TEST2), & + (CCC.GT.TEST1.OR.CCC.GT.TEST2) +#endif + ENDDO + ENDDO + ENDDO + DEALLOCATE(MAXCLA) + ! + !..... counts the fraction of the eliminated configurations + ELIM=(1.D0-DBLE(NCONF)/DBLE(NCONFM))*100.D0 + ! WRITE(994,*) 'NCONF:',NCONF,ELIM + END SUBROUTINE INSNLGQM !/ !/ End of module W3SNL1MD -------------------------------------------- / !/ diff --git a/model/src/w3srcemd.F90 b/model/src/w3srcemd.F90 index 6aa708bb8..a846605d8 100644 --- a/model/src/w3srcemd.F90 +++ b/model/src/w3srcemd.F90 @@ -564,6 +564,7 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & #endif #ifdef W3_NL1 USE W3SNL1MD + USE W3GDATMD, ONLY: IQTPE #endif #ifdef W3_NL2 USE W3SNL2MD @@ -1215,7 +1216,11 @@ SUBROUTINE W3SRCE ( srce_call, IT, ISEA, JSEA, IX, IY, IMOD, & ! 2.b Nonlinear interactions. ! #ifdef W3_NL1 - CALL W3SNL1 ( SPEC, CG1, WNMEAN*DEPTH, VSNL, VDNL ) + IF (IQTPE.GT.0) THEN + CALL W3SNL1 ( SPEC, CG1, WNMEAN*DEPTH, VSNL, VDNL ) + ELSE + CALL W3SNLGQM ( SPEC, CG1, WN1, DEPTH, VSNL, VDNL ) + END IF #endif #ifdef W3_NL2 CALL W3SNL2 ( SPEC, CG1, DEPTH, VSNL, VDNL ) diff --git a/model/src/ww3_ounp.F90 b/model/src/ww3_ounp.F90 index 409888da1..0cbc0b653 100644 --- a/model/src/ww3_ounp.F90 +++ b/model/src/ww3_ounp.F90 @@ -1547,6 +1547,7 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER) #endif #ifdef W3_NL1 USE W3SNL1MD + USE W3GDATMD, ONLY: IQTPE #endif #ifdef W3_NL2 USE W3SNL2MD @@ -2421,7 +2422,11 @@ SUBROUTINE W3EXNC(I,NCID,NREQ,INDREQ,ORDER) END IF IF ( FLSRCE(3) ) THEN #ifdef W3_NL1 - CALL W3SNL1 ( A, CG, WNMEAN*DEPTH, XNL, DIA ) + IF (IQTPE.GT.0) THEN + CALL W3SNL1 ( A, CG, WNMEAN*DEPTH, XNL, DIA ) + ELSE + CALL W3SNLGQM ( A, CG, WN, DEPTH, XNL, DIA ) + END IF #endif #ifdef W3_NL2 CALL W3SNL2 ( A, CG, DEPTH, XNL, DIA ) diff --git a/model/src/ww3_outp.F90 b/model/src/ww3_outp.F90 index a95ec2e93..6d750687a 100644 --- a/model/src/ww3_outp.F90 +++ b/model/src/ww3_outp.F90 @@ -1983,7 +1983,11 @@ SUBROUTINE W3EXPO END IF IF ( FLSRCE(3) ) THEN #ifdef W3_NL1 - CALL W3SNL1 ( A, CG, WNMEAN*DEPTH, XNL, DIA ) + IF (IQTPE.GT.0) THEN + CALL W3SNL1 ( A, CG, WNMEAN*DEPTH, XNL, DIA ) + ELSE + CALL W3SNLGQM ( A, CG, WN, DEPTH, XNL, DIA ) + END IF #endif #ifdef W3_NL2 CALL W3SNL2 ( A, CG, DEPTH, XNL, DIA ) diff --git a/model/src/ww3_trnc.F90 b/model/src/ww3_trnc.F90 index c2049751c..ec69db4dc 100644 --- a/model/src/ww3_trnc.F90 +++ b/model/src/ww3_trnc.F90 @@ -71,6 +71,9 @@ PROGRAM W3TRNC !/ ------------------------------------------------------------------- / USE CONSTANTS +#ifdef W3_NL1 + USE W3ADATMD, ONLY : W3NAUX, W3SETA +#endif USE W3GDATMD, ONLY : W3NMOD, W3SETG, FLAGLL, XFR, GNAME USE W3ODATMD, ONLY : W3NOUT, W3SETO, FNMPRE USE W3SERVMD, ONLY : ITRACE, NEXTLN, EXTCDE @@ -131,6 +134,10 @@ PROGRAM W3TRNC ! CALL W3NMOD ( 1, 6, 6 ) CALL W3SETG ( 1, 6, 6 ) +#ifdef W3_NL1 + CALL W3NAUX ( 6, 6 ) + CALL W3SETA ( 1, 6, 6 ) +#endif CALL W3NOUT ( 6, 6 ) CALL W3SETO ( 1, 6, 6 ) ! diff --git a/regtests/bin/matrix.base b/regtests/bin/matrix.base index 864583358..3d1d84f16 100755 --- a/regtests/bin/matrix.base +++ b/regtests/bin/matrix.base @@ -912,6 +912,8 @@ echo "$rtst -s ST4_TSA -w work_ST4_TSA $ww3 ww3_ts1" >> matrix.body echo "$rtst -s ST6 -w work_ST6 $ww3 ww3_ts1" >> matrix.body echo "$rtst -w work_NL5 -i input_nl5_matrix $ww3 ww3_ts1" >> matrix.body + echo "$rtst -g ST4_T707 -w work_T707GQM -i input_10ms -N $ww3 ww3_ts1" >> matrix.body + echo "$rtst -g ST4_T713 -w work_T713GQM -i input_10ms -N $ww3 ww3_ts1" >> matrix.body fi # fetch limited growth, no switch sharing here diff --git a/regtests/bin/matrix_cmake_datarmor b/regtests/bin/matrix_cmake_datarmor index 4e635f3c2..16c31e47c 100755 --- a/regtests/bin/matrix_cmake_datarmor +++ b/regtests/bin/matrix_cmake_datarmor @@ -129,7 +129,7 @@ main_dir="`cd $main_dir 1>/dev/null 2>&1 && pwd`" export mpi='$MPI_LAUNCH' # Compile option - opt="-f -N -S -T" + opt="-f -N -S -T -o both" # Base run_test command line export rtst="./bin/run_cmake_test $opt" diff --git a/regtests/ww3_ts1/input/namelists_ST4_T475.nml b/regtests/ww3_ts1/input/namelists_ST4_T475.nml new file mode 100644 index 000000000..e104247aa --- /dev/null +++ b/regtests/ww3_ts1/input/namelists_ST4_T475.nml @@ -0,0 +1,7 @@ +&SIN4 BETAMAX = 1.75, SWELLF = 0.66, TAUWSHELTER = 0.3, + SWELLF3 = 0.022, SWELLF4 = 115000.0, SWELLF7 = 432000.00 / +&SDS4 FXFM3 = 2.5 / +&SIC2 IC2ROUGH = 0.001000, IC2VISC = 2.000, IC2DMAX =0.300 / +&SIS2 ISC1 =0.200E+00, IS2BREAK = T, IS2DUPDATE = F, IS2CREEPB = 0.200E+08 / + +END OF NAMELISTS diff --git a/regtests/ww3_ts1/input/namelists_ST4_T700.nml b/regtests/ww3_ts1/input/namelists_ST4_T700.nml index aa6ecdf70..b47cc70d0 100644 --- a/regtests/ww3_ts1/input/namelists_ST4_T700.nml +++ b/regtests/ww3_ts1/input/namelists_ST4_T700.nml @@ -1,4 +1,4 @@ &SDS4 SDSBCHOICE=3, SDSC2 = -3.8, SDSBR = 0.005, - SDSSTRAIN =0., SDSSTRAIN2 = 0., FXFM3 = 20., SDSFACMTF = 400., + FXFM3 = 20., SDSFACMTF = 400., SDSCUM=0., SDSC5 =0. / END OF NAMELISTS diff --git a/regtests/ww3_ts1/input/namelists_ST4_T702.nml b/regtests/ww3_ts1/input/namelists_ST4_T702.nml new file mode 100644 index 000000000..444c02e29 --- /dev/null +++ b/regtests/ww3_ts1/input/namelists_ST4_T702.nml @@ -0,0 +1,12 @@ +&SIN4 BETAMAX = 1.7, SWELLF = 0.60, TAUWSHELTER = 0.2, + SWELLF3 = 0.022, SWELLF4 = 115000.0, SWELLF7 = 432000.00 / +&SDS4 SDSBCHOICE = 3, SDSC2 = -3.80, FXFM3 = 20.00, + SDSBR = 0.005, SDSBT = 0.0011, SDSCUM = 0.300, SDSC5 = 1.0, + SDSMWD = 0.00, SDSFACMTF = 400 / +&SNL1 NLPROP = 25000000.0 / +&SIC2 IC2ROUGH = 0.001000, IC2VISC = 2.000, IC2DMAX =0.300 / +&SIS2 ISC1 =0.200E+00, IS2BREAK = T, IS2DUPDATE = F, IS2CREEPB = 0.200E+08 / +&MISC ICEHINIT = 0.5, ICEHMIN = 0.1, CICE0 = 0.25, NOSW =6, + CICEN = 2.00, LICE = 40000., FLAGTR = 4, FACBERG = 0.2 , + WCOR1=21., WCOR2=0.5 / +END OF NAMELISTS diff --git a/regtests/ww3_ts1/input/namelists_ST4_T707.nml b/regtests/ww3_ts1/input/namelists_ST4_T707.nml new file mode 100644 index 000000000..16f81517d --- /dev/null +++ b/regtests/ww3_ts1/input/namelists_ST4_T707.nml @@ -0,0 +1,13 @@ + &SNL1 IQTYPE = -2, GQMNF1 = 11, GQMNT1 = 6, GQMNQ_OM2 = 6, + TAILNL=-5.0, GQMTHRSAT=5E-5, GQMTHRCOU = 0.05, GQAMP1=1., + GQAMP2=0.0022, GQAMP3=1., GQAMP4=1.0 / + &SIN4 BETAMAX = 1.6, TAUWSHELTER = 0.0 / + &SDS4 SDSBCHOICE=3, SDSC2 = -2.3, SDSBR = 0.005, + FXFM3 = 20, SDSFACMTF = 400., + SDSMWD = 2., SDSCUM = 0.35, SDSNUW =0, SDSC5=1., SDSBRF1=0.5 / +&SIC2 IC2ROUGH = 0.001000, IC2VISC = 2.000, IC2DMAX =0.300 / +&SIS2 ISC1 =0.200E+00, IS2BREAK = T, IS2DUPDATE = F, IS2CREEPB = 0.200E+08 / +&MISC ICEHINIT = 0.5, ICEHMIN = 0.1, CICE0 = 0.25, NOSW =6, + CICEN = 2.00, LICE = 40000., FLAGTR = 4, FACBERG = 0.2 , + WCOR1=21., WCOR2=0.5 / +END OF NAMELISTS diff --git a/regtests/ww3_ts1/input/ww3_grid_ST4_T475.nml b/regtests/ww3_ts1/input/ww3_grid_ST4_T475.nml new file mode 100644 index 000000000..7987e9528 --- /dev/null +++ b/regtests/ww3_ts1/input/ww3_grid_ST4_T475.nml @@ -0,0 +1,225 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III - ww3_grid.nml - Grid pre-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the spectrum parameterization via SPECTRUM_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! SPECTRUM%XFR = 0. ! frequency increment +! SPECTRUM%FREQ1 = 0. ! first frequency (Hz) +! SPECTRUM%NK = 0 ! number of frequencies (wavenumbers) +! SPECTRUM%NTH = 0 ! number of direction bins +! SPECTRUM%THOFF = 0. ! relative offset of first direction [-0.5,0.5] +! -------------------------------------------------------------------- ! +&SPECTRUM_NML + SPECTRUM%XFR = 1.10 + SPECTRUM%FREQ1 = 0.0485 + SPECTRUM%NK = 36 + SPECTRUM%NTH = 24 +/ + +! -------------------------------------------------------------------- ! +! Define the run parameterization via RUN_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! RUN%FLDRY = F ! dry run (I/O only, no calculation) +! RUN%FLCX = F ! x-component of propagation +! RUN%FLCY = F ! y-component of propagation +! RUN%FLCTH = F ! direction shift +! RUN%FLCK = F ! wavenumber shift +! RUN%FLSOU = F ! source terms +! -------------------------------------------------------------------- ! +&RUN_NML + RUN%FLSOU = T +/ + +! -------------------------------------------------------------------- ! +! Define the timesteps parameterization via TIMESTEPS_NML namelist +! +! * It is highly recommended to set up time steps which are multiple +! between them. +! +! * The first time step to calculate is the maximum CFL time step +! which depend on the lowest frequency FREQ1 previously set up and the +! lowest spatial grid resolution in meters DXY. +! reminder : 1 degree=60minutes // 1minute=1mile // 1mile=1.852km +! The formula for the CFL time is : +! Tcfl = DXY / (G / (FREQ1*4*Pi) ) with the constants Pi=3,14 and G=9.8m/s²; +! DTXY ~= 90% Tcfl +! DTMAX ~= 3 * DTXY (maximum global time step limit) +! +! * The refraction time step depends on how strong can be the current velocities +! on your grid : +! DTKTH ~= DTMAX / 2 ! in case of no or light current velocities +! DTKTH ~= DTMAX / 10 ! in case of strong current velocities +! +! * The source terms time step is usually defined between 5s and 60s. +! A common value is 10s. +! DTMIN ~= 10 +! +! * namelist must be terminated with / +! * definitions & defaults: +! TIMESTEPS%DTMAX = 0. ! maximum global time step (s) +! TIMESTEPS%DTXY = 0. ! maximum CFL time step for x-y (s) +! TIMESTEPS%DTKTH = 0. ! maximum CFL time step for k-th (s) +! TIMESTEPS%DTMIN = 0. ! minimum source term time step (s) +! -------------------------------------------------------------------- ! +&TIMESTEPS_NML + TIMESTEPS%DTMAX = 900. + TIMESTEPS%DTXY = 900. + TIMESTEPS%DTKTH = 900. + TIMESTEPS%DTMIN = 15. +/ + +! -------------------------------------------------------------------- ! +! Define the grid to preprocess via GRID_NML namelist +! +! * the tunable parameters for source terms, propagation schemes, and +! numerics are read using namelists. +! * Any namelist found in the folowing sections is temporarily written +! to param.scratch, and read from there if necessary. +! * The order of the namelists is immaterial. +! * Namelists not needed for the given switch settings will be skipped +! automatically +! +! * grid type can be : +! 'RECT' : rectilinear +! 'CURV' : curvilinear +! 'UNST' : unstructured (triangle-based) +! +! * coordinate system can be : +! 'SPHE' : Spherical (degrees) +! 'CART' : Cartesian (meters) +! +! * grid closure can only be applied in spherical coordinates +! +! * grid closure can be : +! 'NONE' : No closure is applied +! 'SMPL' : Simple grid closure. Grid is periodic in the +! : i-index and wraps at i=NX+1. In other words, +! : (NX+1,J) => (1,J). A grid with simple closure +! : may be rectilinear or curvilinear. +! 'TRPL' : Tripole grid closure : Grid is periodic in the +! : i-index and wraps at i=NX+1 and has closure at +! : j=NY+1. In other words, (NX+1,J<=NY) => (1,J) +! : and (I,NY+1) => (NX-I+1,NY). Tripole +! : grid closure requires that NX be even. A grid +! : with tripole closure must be curvilinear. +! +! * The coastline limit depth is the value which distinguish the sea +! points to the land points. All the points with depth values (ZBIN) +! greater than this limit (ZLIM) will be considered as excluded points +! and will never be wet points, even if the water level grows over. +! It can only overwrite the status of a sea point to a land point. +! The value must have a negative value under the mean sea level +! +! * The minimum water depth allowed to compute the model is the absolute +! depth value (DMIN) used in the model if the input depth is lower to +! avoid the model to blow up. +! +! * namelist must be terminated with / +! * definitions & defaults: +! GRID%NAME = 'unset' ! grid name (30 char) +! GRID%NML = 'namelists.nml' ! namelists filename +! GRID%TYPE = 'unset' ! grid type +! GRID%COORD = 'unset' ! coordinate system +! GRID%CLOS = 'unset' ! grid closure +! +! GRID%ZLIM = 0. ! coastline limit depth (m) +! GRID%DMIN = 0. ! abs. minimum water depth (m) +! -------------------------------------------------------------------- ! +&GRID_NML + GRID%NAME = 'HOMOGENEOUS SOURCE TERM TEST' + GRID%NML = '../input/namelists_ST4_T475.nml' + GRID%TYPE = 'RECT' + GRID%COORD = 'SPHE' + GRID%CLOS = 'NONE' + GRID%ZLIM = -5. + GRID%DMIN = 5.75 +/ + +! -------------------------------------------------------------------- ! +! Define the rectilinear grid type via RECT_NML namelist +! - only for RECT grids - +! +! * The minimum grid size is 3x3. +! +! * If the grid increments SX and SY are given in minutes of arc, the scaling +! factor SF must be set to 60. to provide an increment factor in degree. +! +! * If CSTRG='SMPL', then SX is forced to 360/NX. +! +! * value <= value_read / scale_fac +! +! * namelist must be terminated with / +! * definitions & defaults: +! RECT%NX = 0 ! number of points along x-axis +! RECT%NY = 0 ! number of points along y-axis +! +! RECT%SX = 0. ! grid increment along x-axis +! RECT%SY = 0. ! grid increment along y-axis +! RECT%SF = 1. ! scaling division factor for x-y axis +! +! RECT%X0 = 0. ! x-coordinate of lower-left corner (deg) +! RECT%Y0 = 0. ! y-coordinate of lower-left corner (deg) +! RECT%SF0 = 1. ! scaling division factor for x0,y0 coord +! -------------------------------------------------------------------- ! +&RECT_NML + RECT%NX = 3 + RECT%NY = 3 + RECT%SX = 1. + RECT%SY = 1. + RECT%SF = 1.E-2 + RECT%X0 = -1. + RECT%Y0 = -1. + RECT%SF0 = 1.E-2 +/ + +! -------------------------------------------------------------------- ! +! Define the depth to preprocess via DEPTH_NML namelist +! - for RECT and CURV grids - +! +! * if no obstruction subgrid, need to set &MISC FLAGTR = 0 +! +! * The depth value must have negative values under the mean sea level +! +! * value <= value_read * scale_fac +! +! * IDLA : Layout indicator : +! 1 : Read line-by-line bottom to top. (default) +! 2 : Like 1, single read statement. +! 3 : Read line-by-line top to bottom. +! 4 : Like 3, single read statement. +! * IDFM : format indicator : +! 1 : Free format. (default) +! 2 : Fixed format. +! 3 : Unformatted. +! * FORMAT : element format to read : +! '(....)' : auto detected (default) +! '(f10.6)' : float type +! +! * Example : +! IDF SF IDLA IDFM FORMAT FILENAME +! 50 0.001 1 1 '(....)' 'GLOB-30M.bot' +! +! * namelist must be terminated with / +! * definitions & defaults: +! DEPTH%SF = 1. ! scale factor +! DEPTH%FILENAME = 'unset' ! filename +! DEPTH%IDF = 50 ! file unit number +! DEPTH%IDLA = 1 ! layout indicator +! DEPTH%IDFM = 1 ! format indicator +! DEPTH%FORMAT = '(....)' ! formatted read format +! -------------------------------------------------------------------- ! +&DEPTH_NML + DEPTH%SF = -2500. + DEPTH%FILENAME = '../input/HOMOGENEOUS.depth' + DEPTH%IDLA = 3 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input/ww3_grid_ST4_T702.nml b/regtests/ww3_ts1/input/ww3_grid_ST4_T702.nml new file mode 100644 index 000000000..48135e1d9 --- /dev/null +++ b/regtests/ww3_ts1/input/ww3_grid_ST4_T702.nml @@ -0,0 +1,225 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III - ww3_grid.nml - Grid pre-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the spectrum parameterization via SPECTRUM_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! SPECTRUM%XFR = 0. ! frequency increment +! SPECTRUM%FREQ1 = 0. ! first frequency (Hz) +! SPECTRUM%NK = 0 ! number of frequencies (wavenumbers) +! SPECTRUM%NTH = 0 ! number of direction bins +! SPECTRUM%THOFF = 0. ! relative offset of first direction [-0.5,0.5] +! -------------------------------------------------------------------- ! +&SPECTRUM_NML + SPECTRUM%XFR = 1.10 + SPECTRUM%FREQ1 = 0.0485 + SPECTRUM%NK = 36 + SPECTRUM%NTH = 24 +/ + +! -------------------------------------------------------------------- ! +! Define the run parameterization via RUN_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! RUN%FLDRY = F ! dry run (I/O only, no calculation) +! RUN%FLCX = F ! x-component of propagation +! RUN%FLCY = F ! y-component of propagation +! RUN%FLCTH = F ! direction shift +! RUN%FLCK = F ! wavenumber shift +! RUN%FLSOU = F ! source terms +! -------------------------------------------------------------------- ! +&RUN_NML + RUN%FLSOU = T +/ + +! -------------------------------------------------------------------- ! +! Define the timesteps parameterization via TIMESTEPS_NML namelist +! +! * It is highly recommended to set up time steps which are multiple +! between them. +! +! * The first time step to calculate is the maximum CFL time step +! which depend on the lowest frequency FREQ1 previously set up and the +! lowest spatial grid resolution in meters DXY. +! reminder : 1 degree=60minutes // 1minute=1mile // 1mile=1.852km +! The formula for the CFL time is : +! Tcfl = DXY / (G / (FREQ1*4*Pi) ) with the constants Pi=3,14 and G=9.8m/s²; +! DTXY ~= 90% Tcfl +! DTMAX ~= 3 * DTXY (maximum global time step limit) +! +! * The refraction time step depends on how strong can be the current velocities +! on your grid : +! DTKTH ~= DTMAX / 2 ! in case of no or light current velocities +! DTKTH ~= DTMAX / 10 ! in case of strong current velocities +! +! * The source terms time step is usually defined between 5s and 60s. +! A common value is 10s. +! DTMIN ~= 10 +! +! * namelist must be terminated with / +! * definitions & defaults: +! TIMESTEPS%DTMAX = 0. ! maximum global time step (s) +! TIMESTEPS%DTXY = 0. ! maximum CFL time step for x-y (s) +! TIMESTEPS%DTKTH = 0. ! maximum CFL time step for k-th (s) +! TIMESTEPS%DTMIN = 0. ! minimum source term time step (s) +! -------------------------------------------------------------------- ! +&TIMESTEPS_NML + TIMESTEPS%DTMAX = 900. + TIMESTEPS%DTXY = 900. + TIMESTEPS%DTKTH = 900. + TIMESTEPS%DTMIN = 15. +/ + +! -------------------------------------------------------------------- ! +! Define the grid to preprocess via GRID_NML namelist +! +! * the tunable parameters for source terms, propagation schemes, and +! numerics are read using namelists. +! * Any namelist found in the folowing sections is temporarily written +! to param.scratch, and read from there if necessary. +! * The order of the namelists is immaterial. +! * Namelists not needed for the given switch settings will be skipped +! automatically +! +! * grid type can be : +! 'RECT' : rectilinear +! 'CURV' : curvilinear +! 'UNST' : unstructured (triangle-based) +! +! * coordinate system can be : +! 'SPHE' : Spherical (degrees) +! 'CART' : Cartesian (meters) +! +! * grid closure can only be applied in spherical coordinates +! +! * grid closure can be : +! 'NONE' : No closure is applied +! 'SMPL' : Simple grid closure. Grid is periodic in the +! : i-index and wraps at i=NX+1. In other words, +! : (NX+1,J) => (1,J). A grid with simple closure +! : may be rectilinear or curvilinear. +! 'TRPL' : Tripole grid closure : Grid is periodic in the +! : i-index and wraps at i=NX+1 and has closure at +! : j=NY+1. In other words, (NX+1,J<=NY) => (1,J) +! : and (I,NY+1) => (NX-I+1,NY). Tripole +! : grid closure requires that NX be even. A grid +! : with tripole closure must be curvilinear. +! +! * The coastline limit depth is the value which distinguish the sea +! points to the land points. All the points with depth values (ZBIN) +! greater than this limit (ZLIM) will be considered as excluded points +! and will never be wet points, even if the water level grows over. +! It can only overwrite the status of a sea point to a land point. +! The value must have a negative value under the mean sea level +! +! * The minimum water depth allowed to compute the model is the absolute +! depth value (DMIN) used in the model if the input depth is lower to +! avoid the model to blow up. +! +! * namelist must be terminated with / +! * definitions & defaults: +! GRID%NAME = 'unset' ! grid name (30 char) +! GRID%NML = 'namelists.nml' ! namelists filename +! GRID%TYPE = 'unset' ! grid type +! GRID%COORD = 'unset' ! coordinate system +! GRID%CLOS = 'unset' ! grid closure +! +! GRID%ZLIM = 0. ! coastline limit depth (m) +! GRID%DMIN = 0. ! abs. minimum water depth (m) +! -------------------------------------------------------------------- ! +&GRID_NML + GRID%NAME = 'HOMOGENEOUS SOURCE TERM TEST' + GRID%NML = '../input/namelists_ST4_T702.nml' + GRID%TYPE = 'RECT' + GRID%COORD = 'SPHE' + GRID%CLOS = 'NONE' + GRID%ZLIM = -5. + GRID%DMIN = 5.75 +/ + +! -------------------------------------------------------------------- ! +! Define the rectilinear grid type via RECT_NML namelist +! - only for RECT grids - +! +! * The minimum grid size is 3x3. +! +! * If the grid increments SX and SY are given in minutes of arc, the scaling +! factor SF must be set to 60. to provide an increment factor in degree. +! +! * If CSTRG='SMPL', then SX is forced to 360/NX. +! +! * value <= value_read / scale_fac +! +! * namelist must be terminated with / +! * definitions & defaults: +! RECT%NX = 0 ! number of points along x-axis +! RECT%NY = 0 ! number of points along y-axis +! +! RECT%SX = 0. ! grid increment along x-axis +! RECT%SY = 0. ! grid increment along y-axis +! RECT%SF = 1. ! scaling division factor for x-y axis +! +! RECT%X0 = 0. ! x-coordinate of lower-left corner (deg) +! RECT%Y0 = 0. ! y-coordinate of lower-left corner (deg) +! RECT%SF0 = 1. ! scaling division factor for x0,y0 coord +! -------------------------------------------------------------------- ! +&RECT_NML + RECT%NX = 3 + RECT%NY = 3 + RECT%SX = 1. + RECT%SY = 1. + RECT%SF = 1.E-2 + RECT%X0 = -1. + RECT%Y0 = -1. + RECT%SF0 = 1.E-2 +/ + +! -------------------------------------------------------------------- ! +! Define the depth to preprocess via DEPTH_NML namelist +! - for RECT and CURV grids - +! +! * if no obstruction subgrid, need to set &MISC FLAGTR = 0 +! +! * The depth value must have negative values under the mean sea level +! +! * value <= value_read * scale_fac +! +! * IDLA : Layout indicator : +! 1 : Read line-by-line bottom to top. (default) +! 2 : Like 1, single read statement. +! 3 : Read line-by-line top to bottom. +! 4 : Like 3, single read statement. +! * IDFM : format indicator : +! 1 : Free format. (default) +! 2 : Fixed format. +! 3 : Unformatted. +! * FORMAT : element format to read : +! '(....)' : auto detected (default) +! '(f10.6)' : float type +! +! * Example : +! IDF SF IDLA IDFM FORMAT FILENAME +! 50 0.001 1 1 '(....)' 'GLOB-30M.bot' +! +! * namelist must be terminated with / +! * definitions & defaults: +! DEPTH%SF = 1. ! scale factor +! DEPTH%FILENAME = 'unset' ! filename +! DEPTH%IDF = 50 ! file unit number +! DEPTH%IDLA = 1 ! layout indicator +! DEPTH%IDFM = 1 ! format indicator +! DEPTH%FORMAT = '(....)' ! formatted read format +! -------------------------------------------------------------------- ! +&DEPTH_NML + DEPTH%SF = -2500. + DEPTH%FILENAME = '../input/HOMOGENEOUS.depth' + DEPTH%IDLA = 3 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input/ww3_grid_ST4_T707.nml b/regtests/ww3_ts1/input/ww3_grid_ST4_T707.nml new file mode 100644 index 000000000..e6ef84a56 --- /dev/null +++ b/regtests/ww3_ts1/input/ww3_grid_ST4_T707.nml @@ -0,0 +1,225 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III - ww3_grid.nml - Grid pre-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the spectrum parameterization via SPECTRUM_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! SPECTRUM%XFR = 0. ! frequency increment +! SPECTRUM%FREQ1 = 0. ! first frequency (Hz) +! SPECTRUM%NK = 0 ! number of frequencies (wavenumbers) +! SPECTRUM%NTH = 0 ! number of direction bins +! SPECTRUM%THOFF = 0. ! relative offset of first direction [-0.5,0.5] +! -------------------------------------------------------------------- ! +&SPECTRUM_NML + SPECTRUM%XFR = 1.10 + SPECTRUM%FREQ1 = 0.0485 + SPECTRUM%NK = 36 + SPECTRUM%NTH = 24 +/ + +! -------------------------------------------------------------------- ! +! Define the run parameterization via RUN_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! RUN%FLDRY = F ! dry run (I/O only, no calculation) +! RUN%FLCX = F ! x-component of propagation +! RUN%FLCY = F ! y-component of propagation +! RUN%FLCTH = F ! direction shift +! RUN%FLCK = F ! wavenumber shift +! RUN%FLSOU = F ! source terms +! -------------------------------------------------------------------- ! +&RUN_NML + RUN%FLSOU = T +/ + +! -------------------------------------------------------------------- ! +! Define the timesteps parameterization via TIMESTEPS_NML namelist +! +! * It is highly recommended to set up time steps which are multiple +! between them. +! +! * The first time step to calculate is the maximum CFL time step +! which depend on the lowest frequency FREQ1 previously set up and the +! lowest spatial grid resolution in meters DXY. +! reminder : 1 degree=60minutes // 1minute=1mile // 1mile=1.852km +! The formula for the CFL time is : +! Tcfl = DXY / (G / (FREQ1*4*Pi) ) with the constants Pi=3,14 and G=9.8m/s²; +! DTXY ~= 90% Tcfl +! DTMAX ~= 3 * DTXY (maximum global time step limit) +! +! * The refraction time step depends on how strong can be the current velocities +! on your grid : +! DTKTH ~= DTMAX / 2 ! in case of no or light current velocities +! DTKTH ~= DTMAX / 10 ! in case of strong current velocities +! +! * The source terms time step is usually defined between 5s and 60s. +! A common value is 10s. +! DTMIN ~= 10 +! +! * namelist must be terminated with / +! * definitions & defaults: +! TIMESTEPS%DTMAX = 0. ! maximum global time step (s) +! TIMESTEPS%DTXY = 0. ! maximum CFL time step for x-y (s) +! TIMESTEPS%DTKTH = 0. ! maximum CFL time step for k-th (s) +! TIMESTEPS%DTMIN = 0. ! minimum source term time step (s) +! -------------------------------------------------------------------- ! +&TIMESTEPS_NML + TIMESTEPS%DTMAX = 900. + TIMESTEPS%DTXY = 900. + TIMESTEPS%DTKTH = 900. + TIMESTEPS%DTMIN = 15. +/ + +! -------------------------------------------------------------------- ! +! Define the grid to preprocess via GRID_NML namelist +! +! * the tunable parameters for source terms, propagation schemes, and +! numerics are read using namelists. +! * Any namelist found in the folowing sections is temporarily written +! to param.scratch, and read from there if necessary. +! * The order of the namelists is immaterial. +! * Namelists not needed for the given switch settings will be skipped +! automatically +! +! * grid type can be : +! 'RECT' : rectilinear +! 'CURV' : curvilinear +! 'UNST' : unstructured (triangle-based) +! +! * coordinate system can be : +! 'SPHE' : Spherical (degrees) +! 'CART' : Cartesian (meters) +! +! * grid closure can only be applied in spherical coordinates +! +! * grid closure can be : +! 'NONE' : No closure is applied +! 'SMPL' : Simple grid closure. Grid is periodic in the +! : i-index and wraps at i=NX+1. In other words, +! : (NX+1,J) => (1,J). A grid with simple closure +! : may be rectilinear or curvilinear. +! 'TRPL' : Tripole grid closure : Grid is periodic in the +! : i-index and wraps at i=NX+1 and has closure at +! : j=NY+1. In other words, (NX+1,J<=NY) => (1,J) +! : and (I,NY+1) => (NX-I+1,NY). Tripole +! : grid closure requires that NX be even. A grid +! : with tripole closure must be curvilinear. +! +! * The coastline limit depth is the value which distinguish the sea +! points to the land points. All the points with depth values (ZBIN) +! greater than this limit (ZLIM) will be considered as excluded points +! and will never be wet points, even if the water level grows over. +! It can only overwrite the status of a sea point to a land point. +! The value must have a negative value under the mean sea level +! +! * The minimum water depth allowed to compute the model is the absolute +! depth value (DMIN) used in the model if the input depth is lower to +! avoid the model to blow up. +! +! * namelist must be terminated with / +! * definitions & defaults: +! GRID%NAME = 'unset' ! grid name (30 char) +! GRID%NML = 'namelists.nml' ! namelists filename +! GRID%TYPE = 'unset' ! grid type +! GRID%COORD = 'unset' ! coordinate system +! GRID%CLOS = 'unset' ! grid closure +! +! GRID%ZLIM = 0. ! coastline limit depth (m) +! GRID%DMIN = 0. ! abs. minimum water depth (m) +! -------------------------------------------------------------------- ! +&GRID_NML + GRID%NAME = 'HOMOGENEOUS SOURCE TERM TEST' + GRID%NML = '../input/namelists_ST4_T707.nml' + GRID%TYPE = 'RECT' + GRID%COORD = 'SPHE' + GRID%CLOS = 'NONE' + GRID%ZLIM = -5. + GRID%DMIN = 5.75 +/ + +! -------------------------------------------------------------------- ! +! Define the rectilinear grid type via RECT_NML namelist +! - only for RECT grids - +! +! * The minimum grid size is 3x3. +! +! * If the grid increments SX and SY are given in minutes of arc, the scaling +! factor SF must be set to 60. to provide an increment factor in degree. +! +! * If CSTRG='SMPL', then SX is forced to 360/NX. +! +! * value <= value_read / scale_fac +! +! * namelist must be terminated with / +! * definitions & defaults: +! RECT%NX = 0 ! number of points along x-axis +! RECT%NY = 0 ! number of points along y-axis +! +! RECT%SX = 0. ! grid increment along x-axis +! RECT%SY = 0. ! grid increment along y-axis +! RECT%SF = 1. ! scaling division factor for x-y axis +! +! RECT%X0 = 0. ! x-coordinate of lower-left corner (deg) +! RECT%Y0 = 0. ! y-coordinate of lower-left corner (deg) +! RECT%SF0 = 1. ! scaling division factor for x0,y0 coord +! -------------------------------------------------------------------- ! +&RECT_NML + RECT%NX = 3 + RECT%NY = 3 + RECT%SX = 1. + RECT%SY = 1. + RECT%SF = 1.E-2 + RECT%X0 = -1. + RECT%Y0 = -1. + RECT%SF0 = 1.E-2 +/ + +! -------------------------------------------------------------------- ! +! Define the depth to preprocess via DEPTH_NML namelist +! - for RECT and CURV grids - +! +! * if no obstruction subgrid, need to set &MISC FLAGTR = 0 +! +! * The depth value must have negative values under the mean sea level +! +! * value <= value_read * scale_fac +! +! * IDLA : Layout indicator : +! 1 : Read line-by-line bottom to top. (default) +! 2 : Like 1, single read statement. +! 3 : Read line-by-line top to bottom. +! 4 : Like 3, single read statement. +! * IDFM : format indicator : +! 1 : Free format. (default) +! 2 : Fixed format. +! 3 : Unformatted. +! * FORMAT : element format to read : +! '(....)' : auto detected (default) +! '(f10.6)' : float type +! +! * Example : +! IDF SF IDLA IDFM FORMAT FILENAME +! 50 0.001 1 1 '(....)' 'GLOB-30M.bot' +! +! * namelist must be terminated with / +! * definitions & defaults: +! DEPTH%SF = 1. ! scale factor +! DEPTH%FILENAME = 'unset' ! filename +! DEPTH%IDF = 50 ! file unit number +! DEPTH%IDLA = 1 ! layout indicator +! DEPTH%IDFM = 1 ! format indicator +! DEPTH%FORMAT = '(....)' ! formatted read format +! -------------------------------------------------------------------- ! +&DEPTH_NML + DEPTH%SF = -2500. + DEPTH%FILENAME = '../input/HOMOGENEOUS.depth' + DEPTH%IDLA = 3 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input_10ms/namelists_ST4_T707.nml b/regtests/ww3_ts1/input_10ms/namelists_ST4_T707.nml new file mode 100644 index 000000000..0458cd775 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/namelists_ST4_T707.nml @@ -0,0 +1,14 @@ + &SNL1 IQTYPE = -2, GQMNF1 = 11, GQMNT1 = 6, GQMNQ_OM2 = 6, + TAILNL=-5.0, GQMTHRSAT=5E-5, GQMTHRCOU = 0.05, GQAMP1=1., + GQAMP2=0.0022, GQAMP3=1., GQAMP4=1.0 / + &SIN4 BETAMAX = 1.6, TAUWSHELTER = 0.0 / + &SDS4 SDSBCHOICE=3, SDSC2 = -2.3, SDSBR = 0.005, + FXFM3 = 20, SDSFACMTF = 400., + SDSMWD = 2., SDSCUM = 0.35, SDSNUW =0, SDSC5=1., SDSBRF1=0.5 / +&SIC2 IC2ROUGH = 0.001000, IC2VISC = 2.000, IC2DMAX =0.300 / +&SIS2 ISC1 =0.200E+00, IS2BREAK = T, IS2DUPDATE = F, IS2CREEPB = 0.200E+08 / +! DO NOT FORGET TO ADD FLAGTR = 4 for real life runs ... +&MISC ICEHINIT = 0.5, ICEHMIN = 0.1, CICE0 = 0.25, NOSW =6, + CICEN = 2.00, LICE = 40000., FACBERG = 0.2 , + WCOR1=21., WCOR2=0.5 / +END OF NAMELISTS diff --git a/regtests/ww3_ts1/input_10ms/namelists_ST4_T713.nml b/regtests/ww3_ts1/input_10ms/namelists_ST4_T713.nml new file mode 100644 index 000000000..878604430 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/namelists_ST4_T713.nml @@ -0,0 +1,14 @@ + &SNL1 IQTYPE = -2, GQMNF1 = 11, GQMNT1 = 6, GQMNQ_OM2 = 6, + TAILNL=-5.0, GQMTHRSAT=5E-5, GQMTHRCOU = 0.05, GQAMP1=1., + GQAMP2=0.0022, GQAMP3=2. / +&SIN4 BETAMAX = 1.1, TAUWSHELTER = 0.0 / +&SDS4 SDSBCHOICE=3, SDSC2 = -2.5, SDSBR = 0.005, + SDSSTRAIN2 =1.,SDSCUMP=1., FXFM3 = 20, SDSFACMTF = 200., + SDSMWD = 0.9, SDSCUM = 0.3, SDSNUW =0, SDSC5=0.5, SDSBRF1=0.5 / +&SIC2 IC2ROUGH = 0.001000, IC2VISC = 2.000, IC2DMAX =0.300 / +&SIS2 ISC1 =0.200E+00, IS2BREAK = T, IS2DUPDATE = F, IS2CREEPB = 0.200E+08 / +! DO NOT FORGET TO ADD FLAGTR = 4 for real life runs ... +&MISC ICEHINIT = 0.5, ICEHMIN = 0.1, CICE0 = 0.25, NOSW =6, + CICEN = 2.00, LICE = 40000., FACBERG = 0.2 , + WCOR1=21., WCOR2=0.5 / +END OF NAMELISTS diff --git a/regtests/ww3_ts1/input_10ms/points.list b/regtests/ww3_ts1/input_10ms/points.list new file mode 100644 index 000000000..5ad8fde50 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/points.list @@ -0,0 +1 @@ +0.0 0.0 'The_point' diff --git a/regtests/ww3_ts1/input_10ms/switch b/regtests/ww3_ts1/input_10ms/switch new file mode 100644 index 000000000..c3b8938ee --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/switch @@ -0,0 +1 @@ +NOGRB SHRD PR0 FLX0 LN1 ST4 NL1 BT1 DB1 TR0 BS0 IC0 IS0 REF0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 O10 O11 diff --git a/regtests/ww3_ts1/input_10ms/switch_ST4 b/regtests/ww3_ts1/input_10ms/switch_ST4 new file mode 100644 index 000000000..c3b8938ee --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/switch_ST4 @@ -0,0 +1 @@ +NOGRB SHRD PR0 FLX0 LN1 ST4 NL1 BT1 DB1 TR0 BS0 IC0 IS0 REF0 WNT1 WNX1 CRT1 CRX1 O0 O1 O2 O3 O4 O5 O6 O7 O10 O11 diff --git a/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T707.nml b/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T707.nml new file mode 100644 index 000000000..5378ebec3 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T707.nml @@ -0,0 +1,225 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III - ww3_grid.nml - Grid pre-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the spectrum parameterization via SPECTRUM_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! SPECTRUM%XFR = 0. ! frequency increment +! SPECTRUM%FREQ1 = 0. ! first frequency (Hz) +! SPECTRUM%NK = 0 ! number of frequencies (wavenumbers) +! SPECTRUM%NTH = 0 ! number of direction bins +! SPECTRUM%THOFF = 0. ! relative offset of first direction [-0.5,0.5] +! -------------------------------------------------------------------- ! +&SPECTRUM_NML + SPECTRUM%XFR = 1.10 + SPECTRUM%FREQ1 = 0.034 + SPECTRUM%NK = 36 + SPECTRUM%NTH = 36 +/ + +! -------------------------------------------------------------------- ! +! Define the run parameterization via RUN_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! RUN%FLDRY = F ! dry run (I/O only, no calculation) +! RUN%FLCX = F ! x-component of propagation +! RUN%FLCY = F ! y-component of propagation +! RUN%FLCTH = F ! direction shift +! RUN%FLCK = F ! wavenumber shift +! RUN%FLSOU = F ! source terms +! -------------------------------------------------------------------- ! +&RUN_NML + RUN%FLSOU = T +/ + +! -------------------------------------------------------------------- ! +! Define the timesteps parameterization via TIMESTEPS_NML namelist +! +! * It is highly recommended to set up time steps which are multiple +! between them. +! +! * The first time step to calculate is the maximum CFL time step +! which depend on the lowest frequency FREQ1 previously set up and the +! lowest spatial grid resolution in meters DXY. +! reminder : 1 degree=60minutes // 1minute=1mile // 1mile=1.852km +! The formula for the CFL time is : +! Tcfl = DXY / (G / (FREQ1*4*Pi) ) with the constants Pi=3,14 and G=9.8m/s²; +! DTXY ~= 90% Tcfl +! DTMAX ~= 3 * DTXY (maximum global time step limit) +! +! * The refraction time step depends on how strong can be the current velocities +! on your grid : +! DTKTH ~= DTMAX / 2 ! in case of no or light current velocities +! DTKTH ~= DTMAX / 10 ! in case of strong current velocities +! +! * The source terms time step is usually defined between 5s and 60s. +! A common value is 10s. +! DTMIN = 10 +! +! * namelist must be terminated with / +! * definitions & defaults: +! TIMESTEPS%DTMAX = 0. ! maximum global time step (s) +! TIMESTEPS%DTXY = 0. ! maximum CFL time step for x-y (s) +! TIMESTEPS%DTKTH = 0. ! maximum CFL time step for k-th (s) +! TIMESTEPS%DTMIN = 0. ! minimum source term time step (s) +! -------------------------------------------------------------------- ! +&TIMESTEPS_NML + TIMESTEPS%DTMAX = 900. + TIMESTEPS%DTXY = 900. + TIMESTEPS%DTKTH = 900. + TIMESTEPS%DTMIN = 5. +/ + +! -------------------------------------------------------------------- ! +! Define the grid to preprocess via GRID_NML namelist +! +! * the tunable parameters for source terms, propagation schemes, and +! numerics are read using namelists. +! * Any namelist found in the folowing sections is temporarily written +! to param.scratch, and read from there if necessary. +! * The order of the namelists is immaterial. +! * Namelists not needed for the given switch settings will be skipped +! automatically +! +! * grid type can be : +! 'RECT' : rectilinear +! 'CURV' : curvilinear +! 'UNST' : unstructured (triangle-based) +! +! * coordinate system can be : +! 'SPHE' : Spherical (degrees) +! 'CART' : Cartesian (meters) +! +! * grid closure can only be applied in spherical coordinates +! +! * grid closure can be : +! 'NONE' : No closure is applied +! 'SMPL' : Simple grid closure. Grid is periodic in the +! : i-index and wraps at i=NX+1. In other words, +! : (NX+1,J) => (1,J). A grid with simple closure +! : may be rectilinear or curvilinear. +! 'TRPL' : Tripole grid closure : Grid is periodic in the +! : i-index and wraps at i=NX+1 and has closure at +! : j=NY+1. In other words, (NX+1,J<=NY) => (1,J) +! : and (I,NY+1) => (NX-I+1,NY). Tripole +! : grid closure requires that NX be even. A grid +! : with tripole closure must be curvilinear. +! +! * The coastline limit depth is the value which distinguish the sea +! points to the land points. All the points with depth values (ZBIN) +! greater than this limit (ZLIM) will be considered as excluded points +! and will never be wet points, even if the water level grows over. +! It can only overwrite the status of a sea point to a land point. +! The value must have a negative value under the mean sea level +! +! * The minimum water depth allowed to compute the model is the absolute +! depth value (DMIN) used in the model if the input depth is lower to +! avoid the model to blow up. +! +! * namelist must be terminated with / +! * definitions & defaults: +! GRID%NAME = 'unset' ! grid name (30 char) +! GRID%NML = 'namelists.nml' ! namelists filename +! GRID%TYPE = 'unset' ! grid type +! GRID%COORD = 'unset' ! coordinate system +! GRID%CLOS = 'unset' ! grid closure +! +! GRID%ZLIM = 0. ! coastline limit depth (m) +! GRID%DMIN = 0. ! abs. minimum water depth (m) +! -------------------------------------------------------------------- ! +&GRID_NML + GRID%NAME = 'HOMOGENEOUS SOURCE TERM TEST' + GRID%NML = '../input_10ms/namelists_ST4_T707.nml' + GRID%TYPE = 'RECT' + GRID%COORD = 'SPHE' + GRID%CLOS = 'NONE' + GRID%ZLIM = -5. + GRID%DMIN = 5.75 +/ + +! -------------------------------------------------------------------- ! +! Define the rectilinear grid type via RECT_NML namelist +! - only for RECT grids - +! +! * The minimum grid size is 3x3. +! +! * If the grid increments SX and SY are given in minutes of arc, the scaling +! factor SF must be set to 60. to provide an increment factor in degree. +! +! * If CSTRG='SMPL', then SX is forced to 360/NX. +! +! * value <= value_read / scale_fac +! +! * namelist must be terminated with / +! * definitions & defaults: +! RECT%NX = 0 ! number of points along x-axis +! RECT%NY = 0 ! number of points along y-axis +! +! RECT%SX = 0. ! grid increment along x-axis +! RECT%SY = 0. ! grid increment along y-axis +! RECT%SF = 1. ! scaling division factor for x-y axis +! +! RECT%X0 = 0. ! x-coordinate of lower-left corner (deg) +! RECT%Y0 = 0. ! y-coordinate of lower-left corner (deg) +! RECT%SF0 = 1. ! scaling division factor for x0,y0 coord +! -------------------------------------------------------------------- ! +&RECT_NML + RECT%NX = 3 + RECT%NY = 3 + RECT%SX = 1. + RECT%SY = 1. + RECT%SF = 1.E-2 + RECT%X0 = -1. + RECT%Y0 = -1. + RECT%SF0 = 1.E-2 +/ + +! -------------------------------------------------------------------- ! +! Define the depth to preprocess via DEPTH_NML namelist +! - for RECT and CURV grids - +! +! * if no obstruction subgrid, need to set &MISC FLAGTR = 0 +! +! * The depth value must have negative values under the mean sea level +! +! * value <= value_read * scale_fac +! +! * IDLA : Layout indicator : +! 1 : Read line-by-line bottom to top. (default) +! 2 : Like 1, single read statement. +! 3 : Read line-by-line top to bottom. +! 4 : Like 3, single read statement. +! * IDFM : format indicator : +! 1 : Free format. (default) +! 2 : Fixed format. +! 3 : Unformatted. +! * FORMAT : element format to read : +! '(....)' : auto detected (default) +! '(f10.6)' : float type +! +! * Example : +! IDF SF IDLA IDFM FORMAT FILENAME +! 50 0.001 1 1 '(....)' 'GLOB-30M.bot' +! +! * namelist must be terminated with / +! * definitions & defaults: +! DEPTH%SF = 1. ! scale factor +! DEPTH%FILENAME = 'unset' ! filename +! DEPTH%IDF = 50 ! file unit number +! DEPTH%IDLA = 1 ! layout indicator +! DEPTH%IDFM = 1 ! format indicator +! DEPTH%FORMAT = '(....)' ! formatted read format +! -------------------------------------------------------------------- ! +&DEPTH_NML + DEPTH%SF = -2500. + DEPTH%FILENAME = '../input/HOMOGENEOUS.depth' + DEPTH%IDLA = 3 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T713.nml b/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T713.nml new file mode 100644 index 000000000..3efd65adf --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/ww3_grid_ST4_T713.nml @@ -0,0 +1,225 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III - ww3_grid.nml - Grid pre-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the spectrum parameterization via SPECTRUM_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! SPECTRUM%XFR = 0. ! frequency increment +! SPECTRUM%FREQ1 = 0. ! first frequency (Hz) +! SPECTRUM%NK = 0 ! number of frequencies (wavenumbers) +! SPECTRUM%NTH = 0 ! number of direction bins +! SPECTRUM%THOFF = 0. ! relative offset of first direction [-0.5,0.5] +! -------------------------------------------------------------------- ! +&SPECTRUM_NML + SPECTRUM%XFR = 1.10 + SPECTRUM%FREQ1 = 0.034 + SPECTRUM%NK = 36 + SPECTRUM%NTH = 36 +/ + +! -------------------------------------------------------------------- ! +! Define the run parameterization via RUN_NML namelist +! +! * namelist must be terminated with / +! * definitions & defaults: +! RUN%FLDRY = F ! dry run (I/O only, no calculation) +! RUN%FLCX = F ! x-component of propagation +! RUN%FLCY = F ! y-component of propagation +! RUN%FLCTH = F ! direction shift +! RUN%FLCK = F ! wavenumber shift +! RUN%FLSOU = F ! source terms +! -------------------------------------------------------------------- ! +&RUN_NML + RUN%FLSOU = T +/ + +! -------------------------------------------------------------------- ! +! Define the timesteps parameterization via TIMESTEPS_NML namelist +! +! * It is highly recommended to set up time steps which are multiple +! between them. +! +! * The first time step to calculate is the maximum CFL time step +! which depend on the lowest frequency FREQ1 previously set up and the +! lowest spatial grid resolution in meters DXY. +! reminder : 1 degree=60minutes // 1minute=1mile // 1mile=1.852km +! The formula for the CFL time is : +! Tcfl = DXY / (G / (FREQ1*4*Pi) ) with the constants Pi=3,14 and G=9.8m/s²; +! DTXY ~= 90% Tcfl +! DTMAX ~= 3 * DTXY (maximum global time step limit) +! +! * The refraction time step depends on how strong can be the current velocities +! on your grid : +! DTKTH ~= DTMAX / 2 ! in case of no or light current velocities +! DTKTH ~= DTMAX / 10 ! in case of strong current velocities +! +! * The source terms time step is usually defined between 5s and 60s. +! A common value is 10s. +! DTMIN = 10 +! +! * namelist must be terminated with / +! * definitions & defaults: +! TIMESTEPS%DTMAX = 0. ! maximum global time step (s) +! TIMESTEPS%DTXY = 0. ! maximum CFL time step for x-y (s) +! TIMESTEPS%DTKTH = 0. ! maximum CFL time step for k-th (s) +! TIMESTEPS%DTMIN = 0. ! minimum source term time step (s) +! -------------------------------------------------------------------- ! +&TIMESTEPS_NML + TIMESTEPS%DTMAX = 900. + TIMESTEPS%DTXY = 900. + TIMESTEPS%DTKTH = 900. + TIMESTEPS%DTMIN = 5. +/ + +! -------------------------------------------------------------------- ! +! Define the grid to preprocess via GRID_NML namelist +! +! * the tunable parameters for source terms, propagation schemes, and +! numerics are read using namelists. +! * Any namelist found in the folowing sections is temporarily written +! to param.scratch, and read from there if necessary. +! * The order of the namelists is immaterial. +! * Namelists not needed for the given switch settings will be skipped +! automatically +! +! * grid type can be : +! 'RECT' : rectilinear +! 'CURV' : curvilinear +! 'UNST' : unstructured (triangle-based) +! +! * coordinate system can be : +! 'SPHE' : Spherical (degrees) +! 'CART' : Cartesian (meters) +! +! * grid closure can only be applied in spherical coordinates +! +! * grid closure can be : +! 'NONE' : No closure is applied +! 'SMPL' : Simple grid closure. Grid is periodic in the +! : i-index and wraps at i=NX+1. In other words, +! : (NX+1,J) => (1,J). A grid with simple closure +! : may be rectilinear or curvilinear. +! 'TRPL' : Tripole grid closure : Grid is periodic in the +! : i-index and wraps at i=NX+1 and has closure at +! : j=NY+1. In other words, (NX+1,J<=NY) => (1,J) +! : and (I,NY+1) => (NX-I+1,NY). Tripole +! : grid closure requires that NX be even. A grid +! : with tripole closure must be curvilinear. +! +! * The coastline limit depth is the value which distinguish the sea +! points to the land points. All the points with depth values (ZBIN) +! greater than this limit (ZLIM) will be considered as excluded points +! and will never be wet points, even if the water level grows over. +! It can only overwrite the status of a sea point to a land point. +! The value must have a negative value under the mean sea level +! +! * The minimum water depth allowed to compute the model is the absolute +! depth value (DMIN) used in the model if the input depth is lower to +! avoid the model to blow up. +! +! * namelist must be terminated with / +! * definitions & defaults: +! GRID%NAME = 'unset' ! grid name (30 char) +! GRID%NML = 'namelists.nml' ! namelists filename +! GRID%TYPE = 'unset' ! grid type +! GRID%COORD = 'unset' ! coordinate system +! GRID%CLOS = 'unset' ! grid closure +! +! GRID%ZLIM = 0. ! coastline limit depth (m) +! GRID%DMIN = 0. ! abs. minimum water depth (m) +! -------------------------------------------------------------------- ! +&GRID_NML + GRID%NAME = 'HOMOGENEOUS SOURCE TERM TEST' + GRID%NML = '../input_10ms/namelists_ST4_T713.nml' + GRID%TYPE = 'RECT' + GRID%COORD = 'SPHE' + GRID%CLOS = 'NONE' + GRID%ZLIM = -5. + GRID%DMIN = 5.75 +/ + +! -------------------------------------------------------------------- ! +! Define the rectilinear grid type via RECT_NML namelist +! - only for RECT grids - +! +! * The minimum grid size is 3x3. +! +! * If the grid increments SX and SY are given in minutes of arc, the scaling +! factor SF must be set to 60. to provide an increment factor in degree. +! +! * If CSTRG='SMPL', then SX is forced to 360/NX. +! +! * value <= value_read / scale_fac +! +! * namelist must be terminated with / +! * definitions & defaults: +! RECT%NX = 0 ! number of points along x-axis +! RECT%NY = 0 ! number of points along y-axis +! +! RECT%SX = 0. ! grid increment along x-axis +! RECT%SY = 0. ! grid increment along y-axis +! RECT%SF = 1. ! scaling division factor for x-y axis +! +! RECT%X0 = 0. ! x-coordinate of lower-left corner (deg) +! RECT%Y0 = 0. ! y-coordinate of lower-left corner (deg) +! RECT%SF0 = 1. ! scaling division factor for x0,y0 coord +! -------------------------------------------------------------------- ! +&RECT_NML + RECT%NX = 3 + RECT%NY = 3 + RECT%SX = 1. + RECT%SY = 1. + RECT%SF = 1.E-2 + RECT%X0 = -1. + RECT%Y0 = -1. + RECT%SF0 = 1.E-2 +/ + +! -------------------------------------------------------------------- ! +! Define the depth to preprocess via DEPTH_NML namelist +! - for RECT and CURV grids - +! +! * if no obstruction subgrid, need to set &MISC FLAGTR = 0 +! +! * The depth value must have negative values under the mean sea level +! +! * value <= value_read * scale_fac +! +! * IDLA : Layout indicator : +! 1 : Read line-by-line bottom to top. (default) +! 2 : Like 1, single read statement. +! 3 : Read line-by-line top to bottom. +! 4 : Like 3, single read statement. +! * IDFM : format indicator : +! 1 : Free format. (default) +! 2 : Fixed format. +! 3 : Unformatted. +! * FORMAT : element format to read : +! '(....)' : auto detected (default) +! '(f10.6)' : float type +! +! * Example : +! IDF SF IDLA IDFM FORMAT FILENAME +! 50 0.001 1 1 '(....)' 'GLOB-30M.bot' +! +! * namelist must be terminated with / +! * definitions & defaults: +! DEPTH%SF = 1. ! scale factor +! DEPTH%FILENAME = 'unset' ! filename +! DEPTH%IDF = 50 ! file unit number +! DEPTH%IDLA = 1 ! layout indicator +! DEPTH%IDFM = 1 ! format indicator +! DEPTH%FORMAT = '(....)' ! formatted read format +! -------------------------------------------------------------------- ! +&DEPTH_NML + DEPTH%SF = -2500. + DEPTH%FILENAME = '../input/HOMOGENEOUS.depth' + DEPTH%IDLA = 3 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input_10ms/ww3_ounf.nml b/regtests/ww3_ts1/input_10ms/ww3_ounf.nml new file mode 100644 index 000000000..716f1df4b --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/ww3_ounf.nml @@ -0,0 +1,29 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III ww3_ounf.nml - Grid output post-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the output fields to postprocess via FIELD_NML namelist +! -------------------------------------------------------------------- ! +&FIELD_NML + FIELD%TIMESTART = '20000101 000000' + FIELD%TIMESTRIDE = '10' + FIELD%TIMECOUNT = '8000' + FIELD%LIST = 'DPT QP QKK WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA UST' + FIELD%PARTITION = '0 1 2' + FIELD%TYPE = 4 +/ + +! -------------------------------------------------------------------- ! +! Define the content of the output file via FILE_NML namelist +! -------------------------------------------------------------------- ! +&FILE_NML + FILE%IX0 = 2 + FILE%IXN = 2 + FILE%IY0 = 2 + FILE%IYN = 2 +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input_10ms/ww3_ounp_spec.nml b/regtests/ww3_ts1/input_10ms/ww3_ounp_spec.nml new file mode 100644 index 000000000..34bac9764 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/ww3_ounp_spec.nml @@ -0,0 +1,48 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III ww3_ounp.nml - Point output post-processing ! +! -------------------------------------------------------------------- ! + +! -------------------------------------------------------------------- ! +! Define the output fields to postprocess via POINT_NML namelist +! -------------------------------------------------------------------- ! +&POINT_NML + POINT%TIMESTART = '20000101 000000' + POINT%TIMESTRIDE = '1800.' + POINT%TIMECOUNT = '1000' + POINT%TIMESPLIT = 4 + POINT%BUFFER = 100 + POINT%TYPE = 3 +/ + +! -------------------------------------------------------------------- ! +! Define the content of the output file via FILE_NML namelist +! -------------------------------------------------------------------- ! +&FILE_NML +/ + +! -------------------------------------------------------------------- ! +! Define the type 0, inventory of file +! -------------------------------------------------------------------- ! + + +! -------------------------------------------------------------------- ! +! Define the type 1, spectra via SPECTRA_NML namelist +! -------------------------------------------------------------------- ! +&SPECTRA_NML +/ + +! -------------------------------------------------------------------- ! +! Define the type 2, mean parameter via PARAM_NML namelist +! -------------------------------------------------------------------- ! +&PARAM_NML +/ + +! -------------------------------------------------------------------- ! +! Define the type 3, source terms via SOURCE_NML namelist +! -------------------------------------------------------------------- ! +&SOURCE_NML +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- ! diff --git a/regtests/ww3_ts1/input_10ms/ww3_shel.nml b/regtests/ww3_ts1/input_10ms/ww3_shel.nml new file mode 100644 index 000000000..d83106ed1 --- /dev/null +++ b/regtests/ww3_ts1/input_10ms/ww3_shel.nml @@ -0,0 +1,54 @@ +! -------------------------------------------------------------------- ! +! WAVEWATCH III ww3_shel.nml - single-grid model ! +! -------------------------------------------------------------------- ! + + +! -------------------------------------------------------------------- ! +! Define top-level model parameters via DOMAIN_NML namelist +! -------------------------------------------------------------------- ! +&DOMAIN_NML + DOMAIN%START = '20000101 000000' + DOMAIN%STOP = '20000106 000000' +/ + +! -------------------------------------------------------------------- ! +! Define each forcing via the INPUT_NML namelist +! -------------------------------------------------------------------- ! +&INPUT_NML + INPUT%FORCING%WINDS = 'H' +/ + +! -------------------------------------------------------------------- ! +! Define the output types point parameters via OUTPUT_TYPE_NML namelist +! -------------------------------------------------------------------- ! +&OUTPUT_TYPE_NML + TYPE%FIELD%LIST = 'DPT QP QKK WND ICE HS MSS MSD FAW WCC WCF WCH WCM FOC TAW CHA UST' + TYPE%POINT%FILE = '../input_10ms/points.list' +/ + +! -------------------------------------------------------------------- ! +! Define output dates via OUTPUT_DATE_NML namelist +! -------------------------------------------------------------------- ! +&OUTPUT_DATE_NML + DATE%FIELD = '19680606 000000' '1800' '20230618 000000' + DATE%POINT = '19680606 000000' '1800' '20230618 000000' +/ + +! -------------------------------------------------------------------- ! +! Define homogeneous input via HOMOG_COUNT_NML and HOMOG_INPUT_NML namelist +! -------------------------------------------------------------------- ! +&HOMOG_COUNT_NML + HOMOG_COUNT%N_CUR = 0 + HOMOG_COUNT%N_WND = 1 +/ + +&HOMOG_INPUT_NML + HOMOG_INPUT(1)%NAME = 'WND' + HOMOG_INPUT(1)%VALUE1 = 10. + HOMOG_INPUT(1)%VALUE2 = 270. + HOMOG_INPUT(1)%VALUE3 = 0. +/ + +! -------------------------------------------------------------------- ! +! WAVEWATCH III - end of namelist ! +! -------------------------------------------------------------------- !