New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11597 for NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex – NEMO

Ignore:
Timestamp:
2019-09-25T20:20:19+02:00 (5 years ago)
Author:
nicolasmartin
Message:

Continuation of coding rules application
Recovery of some sections deleted by the previous commit

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex

    r11596 r11597  
    5353MISC correspond to "extracting tendency terms" or "vorticity balance"?} 
    5454 
     55%% ================================================================================================= 
    5556\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 
    5657\label{sec:DYN_divcur_wzv} 
    5758 
    58 %-------------------------------------------------------------------------------------------------------------- 
    59 %           Horizontal divergence and relative vorticity 
    60 %-------------------------------------------------------------------------------------------------------------- 
     59%% ================================================================================================= 
    6160\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
    6261\label{subsec:DYN_divcur} 
     
    9291the nonlinear advection and of the vertical velocity respectively. 
    9392 
    94 %-------------------------------------------------------------------------------------------------------------- 
    95 %           Sea Surface Height evolution 
    96 %-------------------------------------------------------------------------------------------------------------- 
     93%% ================================================================================================= 
    9794\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
    9895\label{subsec:DYN_sshwzv} 
     
    152149(see \autoref{subsec:DOM_Num_Index_vertical}). 
    153150 
     151 
     152%% ================================================================================================= 
    154153\section{Coriolis and advection: vector invariant form} 
    155154\label{sec:DYN_adv_cor_vect} 
    156 %-----------------------------------------nam_dynadv---------------------------------------------------- 
    157155 
    158156\begin{listing} 
     
    161159  \label{lst:namdyn_adv} 
    162160\end{listing} 
    163 %------------------------------------------------------------------------------------------------------------- 
    164161 
    165162The vector invariant form of the momentum equations is the one most often used in 
     
    172169\autoref{chap:LBC}. 
    173170 
     171%% ================================================================================================= 
    174172\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 
    175173\label{subsec:DYN_vor} 
    176 %------------------------------------------nam_dynvor---------------------------------------------------- 
    177174 
    178175\begin{listing} 
     
    181178  \label{lst:namdyn_vor} 
    182179\end{listing} 
    183 %------------------------------------------------------------------------------------------------------------- 
    184180 
    185181Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables. 
     
    195191The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. 
    196192 
    197 %------------------------------------------------------------- 
    198193%                 enstrophy conserving scheme 
    199 %------------------------------------------------------------- 
     194%% ================================================================================================= 
    200195\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})} 
    201196\label{subsec:DYN_vor_ens} 
     
    218213\end{equation} 
    219214 
    220 %------------------------------------------------------------- 
    221215%                 energy conserving scheme 
    222 %------------------------------------------------------------- 
     216%% ================================================================================================= 
    223217\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})} 
    224218\label{subsec:DYN_vor_ene} 
     
    238232\end{equation} 
    239233 
    240 %------------------------------------------------------------- 
    241234%                 mix energy/enstrophy conserving scheme 
    242 %------------------------------------------------------------- 
     235%% ================================================================================================= 
    243236\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})} 
    244237\label{subsec:DYN_vor_mix} 
     
    263256\] 
    264257 
    265 %------------------------------------------------------------- 
    266258%                 energy and enstrophy conserving scheme 
    267 %------------------------------------------------------------- 
     259%% ================================================================================================= 
    268260\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})} 
    269261\label{subsec:DYN_vor_een} 
     
    353345leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. 
    354346 
    355 %-------------------------------------------------------------------------------------------------------------- 
    356 %           Kinetic Energy Gradient term 
    357 %-------------------------------------------------------------------------------------------------------------- 
     347%% ================================================================================================= 
    358348\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
    359349\label{subsec:DYN_keg} 
     
    373363\] 
    374364 
    375 %-------------------------------------------------------------------------------------------------------------- 
    376 %           Vertical advection term 
    377 %-------------------------------------------------------------------------------------------------------------- 
     365%% ================================================================================================= 
    378366\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} 
    379367\label{subsec:DYN_zad} 
     
    400388an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 
    401389 
     390 
     391%% ================================================================================================= 
    402392\section{Coriolis and advection: flux form} 
    403393\label{sec:DYN_adv_cor_flux} 
    404 %------------------------------------------nam_dynadv---------------------------------------------------- 
    405  
    406 %------------------------------------------------------------------------------------------------------------- 
    407394 
    408395Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables. 
     
    413400no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. 
    414401 
    415 %-------------------------------------------------------------------------------------------------------------- 
    416 %           Coriolis plus curvature metric terms 
    417 %-------------------------------------------------------------------------------------------------------------- 
     402%% ================================================================================================= 
    418403\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 
    419404\label{subsec:DYN_cor_flux} 
     
    434419This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). 
    435420 
    436 %-------------------------------------------------------------------------------------------------------------- 
    437 %           Flux form Advection term 
    438 %-------------------------------------------------------------------------------------------------------------- 
     421%% ================================================================================================= 
    439422\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} 
    440423\label{subsec:DYN_adv_flux} 
     
    467450and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. 
    468451 
    469 %------------------------------------------------------------- 
    470452%                 2nd order centred scheme 
    471 %------------------------------------------------------------- 
     453%% ================================================================================================= 
    472454\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})} 
    473455\label{subsec:DYN_adv_cen2} 
     
    490472so $u$ and $v$ are the \emph{now} velocities. 
    491473 
    492 %------------------------------------------------------------- 
    493474%                 UBS scheme 
    494 %------------------------------------------------------------- 
     475%% ================================================================================================= 
    495476\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})} 
    496477\label{subsec:DYN_adv_ubs} 
     
    543524%%% 
    544525 
     526%% ================================================================================================= 
    545527\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
    546528\label{sec:DYN_hpg} 
    547 %------------------------------------------nam_dynhpg--------------------------------------------------- 
    548529 
    549530\begin{listing} 
     
    552533  \label{lst:namdyn_hpg} 
    553534\end{listing} 
    554 %------------------------------------------------------------------------------------------------------------- 
    555535 
    556536Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables. 
     
    566546At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied. 
    567547 
    568 %-------------------------------------------------------------------------------------------------------------- 
    569 %           z-coordinate with full step 
    570 %-------------------------------------------------------------------------------------------------------------- 
     548%% ================================================================================================= 
    571549\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})} 
    572550\label{subsec:DYN_hpg_zco} 
     
    611589\autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$. 
    612590 
    613 %-------------------------------------------------------------------------------------------------------------- 
    614 %           z-coordinate with partial step 
    615 %-------------------------------------------------------------------------------------------------------------- 
     591%% ================================================================================================= 
    616592\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})} 
    617593\label{subsec:DYN_hpg_zps} 
     
    632608module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}. 
    633609 
    634 %-------------------------------------------------------------------------------------------------------------- 
    635 %           s- and s-z-coordinates 
    636 %-------------------------------------------------------------------------------------------------------------- 
     610%% ================================================================================================= 
    637611\subsection{$S$- and $Z$-$S$-coordinates} 
    638612\label{subsec:DYN_hpg_sco} 
     
    681655This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. 
    682656 
     657%% ================================================================================================= 
    683658\subsection{Ice shelf cavity} 
    684659\label{subsec:DYN_hpg_isf} 
     
    697672\autoref{subsec:DYN_hpg_sco}. 
    698673 
    699 %-------------------------------------------------------------------------------------------------------------- 
    700 %           Time-scheme 
    701 %-------------------------------------------------------------------------------------------------------------- 
     674%% ================================================================================================= 
    702675\subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})} 
    703676\label{subsec:DYN_hpg_imp} 
     
    758731This option is controlled by  \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter. 
    759732 
     733%% ================================================================================================= 
    760734\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})} 
    761735\label{sec:DYN_spg} 
    762 %-----------------------------------------nam_dynspg---------------------------------------------------- 
    763736 
    764737\begin{listing} 
     
    767740  \label{lst:namdyn_spg} 
    768741\end{listing} 
    769 %------------------------------------------------------------------------------------------------------------ 
    770742 
    771743Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables. 
     
    795767As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}. 
    796768 
    797 %-------------------------------------------------------------------------------------------------------------- 
    798 % Explicit free surface formulation 
    799 %-------------------------------------------------------------------------------------------------------------- 
     769%% ================================================================================================= 
    800770\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})} 
    801771\label{subsec:DYN_spg_exp} 
     
    821791Thus, nothing is done in the \mdl{dynspg\_exp} module. 
    822792 
    823 %-------------------------------------------------------------------------------------------------------------- 
    824 % Split-explict free surface formulation 
    825 %-------------------------------------------------------------------------------------------------------------- 
     793%% ================================================================================================= 
    826794\subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})} 
    827795\label{subsec:DYN_spg_ts} 
    828 %------------------------------------------namsplit----------------------------------------------------------- 
    829796% 
    830797%\nlst{namsplit} 
    831 %------------------------------------------------------------------------------------------------------------- 
    832798 
    833799The split-explicit free surface formulation used in \NEMO\ (\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true), 
     
    10651031%>>>>>=============== 
    10661032 
    1067 %-------------------------------------------------------------------------------------------------------------- 
    1068 % Filtered free surface formulation 
    1069 %-------------------------------------------------------------------------------------------------------------- 
     1033%% ================================================================================================= 
    10701034\subsection{Filtered free surface (\forcode{dynspg_flt?})} 
    10711035\label{subsec:DYN_spg_fltp} 
     
    10931057It is computed once and for all and applies to all ocean time steps. 
    10941058 
     1059%% ================================================================================================= 
    10951060\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})} 
    10961061\label{sec:DYN_ldf} 
    1097 %------------------------------------------nam_dynldf---------------------------------------------------- 
    10981062 
    10991063\begin{listing} 
     
    11021066  \label{lst:namdyn_ldf} 
    11031067\end{listing} 
    1104 %------------------------------------------------------------------------------------------------------------- 
    11051068 
    11061069Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables. 
     
    11301093} 
    11311094 
    1132 %           Vertical diffusion term 
    1133 % External Forcing 
    1134 % Wetting and drying 
    1135 % Time evolution term 
     1095%% ================================================================================================= 
     1096\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})} 
     1097\label{subsec:DYN_ldf_lap} 
     1098 
     1099For lateral iso-level diffusion, the discrete operator is: 
     1100\begin{equation} 
     1101  \label{eq:DYN_ldf_lap} 
     1102  \left\{ 
     1103    \begin{aligned} 
     1104      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm} 
     1105          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[ 
     1106        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\ 
     1107      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm} 
     1108          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[ 
     1109        {A_f^{lm} \;e_{3f} \zeta } \right] 
     1110    \end{aligned} 
     1111  \right. 
     1112\end{equation} 
     1113 
     1114As explained in \autoref{subsec:MB_ldf}, 
     1115this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and 
     1116ensures a complete separation between the vorticity and divergence parts of the momentum diffusion. 
     1117 
     1118%% ================================================================================================= 
     1119\subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})} 
     1120\label{subsec:DYN_ldf_iso} 
     1121 
     1122A rotation of the lateral momentum diffusion operator is needed in several cases: 
     1123for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and 
     1124for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or 
     1125geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate. 
     1126In the partial step case, coordinates are horizontal except at the deepest level and 
     1127no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}. 
     1128The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on 
     1129each momentum component. 
     1130It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry. 
     1131The resulting discrete representation is: 
     1132\begin{equation} 
     1133  \label{eq:DYN_ldf_iso} 
     1134  \begin{split} 
     1135    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\ 
     1136    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_T^{lm}  \left( 
     1137              {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u] 
     1138                -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}} 
     1139            \right)} \right]}    \right. \\ 
     1140    & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} 
     1141            }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f} 
     1142            \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}} 
     1143        \right)} \right] \\ 
     1144    &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline 
     1145              {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } 
     1146        \right.} \right. \\ 
     1147    &  \ \qquad \qquad \qquad \quad\ 
     1148    - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\ 
     1149    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 
     1150                +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} 
     1151                \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\ 
     1152    D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\ 
     1153    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_f^{lm}  \left( 
     1154              {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v] 
     1155                -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}} 
     1156            \right)} \right]}    \right. \\ 
     1157    & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 
     1158            }\,\delta_{j} [v] - e_{1t}\, r_{2t} 
     1159            \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}} 
     1160        \right)} \right] \\ 
     1161    & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 
     1162              {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\ 
     1163    &  \ \qquad \qquad \qquad \quad\ 
     1164    - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\ 
     1165    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 
     1166                +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} 
     1167                \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 
     1168  \end{split} 
     1169\end{equation} 
     1170where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and 
     1171the surface of computation ($z$- or $s$-surfaces). 
     1172The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}). 
     1173 
     1174%% ================================================================================================= 
     1175\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})} 
     1176\label{subsec:DYN_ldf_bilap} 
     1177 
     1178The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice. 
     1179It requires an additional assumption on boundary conditions: 
     1180the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, 
     1181while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). 
     1182%%% 
     1183\gmcomment{add a remark on the the change in the position of the coefficient} 
     1184%%% 
     1185 
     1186%% ================================================================================================= 
     1187\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})} 
     1188\label{sec:DYN_zdf} 
     1189 
     1190Options are defined through the \nam{zdf}{zdf} namelist variables. 
     1191The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step. 
     1192Two time stepping schemes can be used for the vertical diffusion term: 
     1193$(a)$ a forward time differencing scheme 
     1194(\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or 
     1195$(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp}) 
     1196(see \autoref{chap:TD}). 
     1197Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics. 
     1198 
     1199The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is. 
     1200The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form: 
     1201\[ 
     1202  % \label{eq:DYN_zdf} 
     1203  \left\{ 
     1204    \begin{aligned} 
     1205      D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} } 
     1206        \ \delta_{k+1/2} [\,u\,]         \right]     \\ 
     1207      \\ 
     1208      D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} } 
     1209        \ \delta_{k+1/2} [\,v\,]         \right] 
     1210    \end{aligned} 
     1211  \right. 
     1212\] 
     1213where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients. 
     1214The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}). 
     1215 
     1216The surface boundary condition on momentum is the stress exerted by the wind. 
     1217At the surface, the momentum fluxes are prescribed as the boundary condition on 
     1218the vertical turbulent momentum fluxes, 
     1219\begin{equation} 
     1220  \label{eq:DYN_zdf_sbc} 
     1221  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1} 
     1222  = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v } 
     1223\end{equation} 
     1224where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in 
     1225the (\textbf{i},\textbf{j}) coordinate system. 
     1226The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in 
     1227the vertical over the mixed layer depth. 
     1228If the vertical mixing coefficient is small (when no mixed layer scheme is used) 
     1229the surface stress enters only the top model level, as a body force. 
     1230The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}). 
     1231 
     1232The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation 
     1233(see \autoref{sec:ZDF_drg}) 
     1234 
     1235%% ================================================================================================= 
     1236\section{External forcings} 
     1237\label{sec:DYN_forcing} 
     1238 
     1239Besides the surface and bottom stresses (see the above section) 
     1240which are introduced as boundary conditions on the vertical mixing, 
     1241three other forcings may enter the dynamical equations by affecting the surface pressure gradient. 
     1242 
     1243(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}), 
     1244the atmospheric pressure is taken into account when computing the surface pressure gradient. 
     1245 
     1246(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}), 
     1247the tidal potential is taken into account when computing the surface pressure gradient. 
     1248 
     1249(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used 
     1250(\ie\ when the sea-ice is embedded in the ocean), 
     1251the snow-ice mass is taken into account when computing the surface pressure gradient. 
     1252 
     1253 
     1254\gmcomment{ missing : the lateral boundary condition !!!   another external forcing 
     1255 } 
     1256 
     1257%% ================================================================================================= 
     1258\section{Wetting and drying } 
     1259\label{sec:DYN_wetdry} 
     1260 
     1261There are two main options for wetting and drying code (wd): 
     1262(a) an iterative limiter (il) and (b) a directional limiter (dl). 
     1263The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO 
     1264MS 
     1265which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative 
     1266limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$ 
     1267and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated 
     1268by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$. 
     1269 
     1270\begin{listing} 
     1271  \nlst{namwad} 
     1272  \caption{\forcode{&namwad}} 
     1273  \label{lst:namwad} 
     1274\end{listing} 
     1275 
     1276The following terminology is used. The depth of the topography (positive downwards) 
     1277at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code. 
     1278The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign 
     1279conventions used, the water depth, $h$, is the height of the free surface plus the depth of the 
     1280topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$). 
     1281 
     1282Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be 
     1283covered by water. They require the topography specified with a model 
     1284configuration to have negative depths at points where the land is higher than the 
     1285topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an 
     1286initial state with zero sea surface height elevation. 
     1287The user can choose to compute the vertical grid and heights in the model relative to 
     1288a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths 
     1289(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays). 
     1290 
     1291Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''. 
     1292$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies 
     1293with very steep slopes require larger values for normal choices of time-step. Surface fluxes 
     1294are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers. 
     1295The fluxes are tappered down using a $\mathrm{tanh}$ weighting function 
     1296to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow. 
     1297The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$. 
     1298The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$. 
     1299 
     1300Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration 
     1301and in ``realistic'' configurations covering parts of the north-west European shelf. 
     1302All these configurations have used pure sigma coordinates. It is expected that 
     1303the wetting and drying code will work in domains with more general s-coordinates provided 
     1304the coordinates are pure sigma in the region where wetting and drying actually occurs. 
     1305 
     1306The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter. 
     1307The final sub-section covers some additional considerations that are relevant to both schemes. 
     1308 
     1309 
     1310%   Iterative limiters 
     1311%% ================================================================================================= 
     1312\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})} 
     1313\label{subsec:DYN_wd_directional_limiter} 
     1314 
     1315The principal idea of the directional limiter is that 
     1316water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn_wdmin1}{rn\_wdmin1}). 
     1317 
     1318All the changes associated with this option are made to the barotropic solver for the non-linear 
     1319free surface code within dynspg\_ts. 
     1320On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells 
     1321and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells 
     1322whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells 
     1323to be at least the minimum depth and hence is able to conserve mass / volume. 
     1324 
     1325The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj). 
     1326If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the 
     1327flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets 
     1328\np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases 
     1329from 2 * \np{rn_wdmin1}{rn\_wdmin1} to \np{rn_wdmin1}{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases. 
     1330 
     1331At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen 
     1332also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have 
     1333chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells. 
     1334The rationale for setting the velocity to zero is that it is the momentum equations that are being solved 
     1335and the total momentum of the upstream cell (treating it as a finite volume) should be considered 
     1336to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its 
     1337treatment in the calculation of the flux of mass across the cell face. 
     1338 
     1339 
     1340\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic 
     1341timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than 
     1342or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer 
     1343fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because 
     1344the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts 
     1345to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the 
     1346baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask. 
     1347 
     1348%   Iterative limiters 
     1349 
     1350%% ================================================================================================= 
     1351\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})} 
     1352\label{subsec:DYN_wd_iterative_limiter} 
     1353 
     1354%% ================================================================================================= 
     1355\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})} 
     1356\label{subsec:DYN_wd_il_spg_limiter} 
     1357 
     1358The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry'' 
     1359or may become dry within the next time-step using an iterative method. 
     1360 
     1361The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows: 
     1362 
     1363The continuity equation for the total water depth in a column 
     1364\begin{equation} 
     1365  \label{eq:DYN_wd_continuity} 
     1366  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 . 
     1367\end{equation} 
     1368can be written in discrete form  as 
     1369 
     1370\begin{align} 
     1371  \label{eq:DYN_wd_continuity_2} 
     1372  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) ) 
     1373  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\ 
     1374  &= \mathrm{zzflx}_{i,j} . 
     1375\end{align} 
     1376 
     1377In the above $h$ is the depth of the water in the column at point $(i,j)$, 
     1378$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and 
     1379$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is 
     1380the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t = 
     1381t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and 
     1382$\mathrm{zzflx}$ is the sum of the fluxes through all the faces. 
     1383 
     1384The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell 
     1385(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly 
     1386 
     1387\begin{equation} 
     1388  \label{eq:DYN_wd_zzflx_p_n_1} 
     1389  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} . 
     1390\end{equation} 
     1391 
     1392The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until 
     1393none of the cells will ``dry out''. To be precise the fluxes are limited until none of the 
     1394cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$. 
     1395 
     1396Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and 
     1397$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients, 
     1398$\mathrm{zcoef}_{i,j}^{(m)}$ such that: 
     1399 
     1400\begin{equation} 
     1401  \label{eq:DYN_wd_continuity_coef} 
     1402  \begin{split} 
     1403    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\ 
     1404    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j} 
     1405  \end{split} 
     1406\end{equation} 
     1407 
     1408where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around 
     1409cells that would otherwise dry. 
     1410 
     1411The iteration is initialised by setting 
     1412 
     1413\begin{equation} 
     1414  \label{eq:DYN_wd_zzflx_initial} 
     1415  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} . 
     1416\end{equation} 
     1417 
     1418The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the 
     1419cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell 
     1420times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this 
     1421condition is 
     1422 
     1423\begin{equation} 
     1424  \label{eq:DYN_wd_continuity_if} 
     1425  h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) . 
     1426\end{equation} 
     1427 
     1428Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum 
     1429outward flux that can be allowed and still maintain the minimum wet depth: 
     1430 
     1431\begin{equation} 
     1432  \label{eq:DYN_wd_max_flux} 
     1433  \begin{split} 
     1434    \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
     1435    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] 
     1436  \end{split} 
     1437\end{equation} 
     1438 
     1439Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is 
     1440this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an 
     1441expression for the coefficient needed to multiply the outward flux at this cell in order 
     1442to avoid drying. 
     1443 
     1444\begin{equation} 
     1445  \label{eq:DYN_wd_continuity_nxtcoef} 
     1446  \begin{split} 
     1447    \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\ 
     1448    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} } 
     1449  \end{split} 
     1450\end{equation} 
     1451 
     1452Only the outward flux components are altered but, of course, outward fluxes from one cell 
     1453are inward fluxes to adjacent cells and the balance in these cells may need subsequent 
     1454adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux 
     1455across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration 
     1456if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then 
     1457the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by 
     1458the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across 
     1459the faces of the cells do not ``compete'' (they do not over-write each other) and one 
     1460would expect the scheme to converge relatively quickly. The scheme is flux based so 
     1461conserves mass. It also conserves constant tracers for the same reason that the 
     1462directional limiter does. 
     1463 
     1464 
     1465%      Surface pressure gradients 
     1466%% ================================================================================================= 
     1467\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})} 
     1468\label{subsec:DYN_wd_il_spg} 
     1469 
     1470At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the 
     1471topography is sloping at these points the sea-surface will have a similar slope and there 
     1472will hence be very large horizontal pressure gradients at these points. The WAD modifies 
     1473the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such 
     1474points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie 
     1475between $0$ and $1$. 
     1476 
     1477We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated 
     1478at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the 
     1479neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals 
     1480variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid 
     1481column.  The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}. 
     1482 
     1483%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     1484\begin{figure}[!ht] 
     1485  \centering 
     1486  \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg} 
     1487  \caption[Combinations controlling the limiting of the horizontal pressure gradient in 
     1488  wetting and drying regimes]{ 
     1489    Three possible combinations of the logical variables controlling the 
     1490    limiting of the horizontal pressure gradient in wetting and drying regimes} 
     1491  \label{fig:DYN_WAD_dynhpg} 
     1492\end{figure} 
     1493%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     1494 
     1495The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at 
     1496both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and 
     1497the minimum height of the sea surface at the two points is greater than the maximum height 
     1498of the topography at the two points: 
     1499 
     1500\begin{equation} 
     1501  \label{eq:DYN_ll_tmp1} 
     1502  \begin{split} 
     1503    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\ 
     1504                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\ 
     1505                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\ 
     1506                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\ 
     1507                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 } 
     1508  \end{split} 
     1509\end{equation} 
     1510 
     1511The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height 
     1512of the sea surface at the two points is greater than the maximum height of the topography 
     1513at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ 
     1514 
     1515\begin{equation} 
     1516  \label{eq:DYN_ll_tmp2} 
     1517  \begin{split} 
     1518    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\ 
     1519    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\ 
     1520    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) . 
     1521  \end{split} 
     1522\end{equation} 
     1523 
     1524If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$ 
     1525point is unmodified. If both logicals are false zhpi is set to zero. 
     1526 
     1527If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure 
     1528gradient is multiplied through by zcpx which is the absolute value of the difference in 
     1529the water depths at the two points divided by the difference in the surface heights at the 
     1530two points. Thus the sign of the sea surface height gradient is retained but the magnitude 
     1531of the pressure force is determined by the difference in water depths rather than the 
     1532difference in surface height between the two points. Note that dividing by the difference 
     1533between the sea surface heights can be problematic if the heights approach parity. An 
     1534additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such 
     1535conditions. 
     1536 
     1537%% ================================================================================================= 
     1538\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})} 
     1539\label{subsec:DYN_WAD_additional} 
     1540 
     1541In the very shallow water where wetting and drying occurs the parametrisation of 
     1542bottom drag is clearly very important. In order to promote stability 
     1543it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach. 
     1544 
     1545Suitable specifcation of the surface heat flux in wetting and drying domains in forced and 
     1546coupled simulations needs further consideration. In order to prevent freezing or boiling 
     1547in uncoupled integrations the net surface heat fluxes need to be appropriately limited. 
     1548 
     1549%      The WAD test cases 
     1550%% ================================================================================================= 
     1551\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})} 
     1552\label{subsec:DYN_WAD_test_cases} 
     1553 
     1554See the WAD tests MY\_DOC documention for details of the WAD test cases. 
     1555 
     1556 
     1557 
     1558%% ================================================================================================= 
     1559\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})} 
     1560\label{sec:DYN_nxt} 
     1561 
     1562 
     1563Options are defined through the \nam{dom}{dom} namelist variables. 
     1564The general framework for dynamics time stepping is a leap-frog scheme, 
     1565\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}). 
     1566The scheme is applied to the velocity, except when 
     1567using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux}) 
     1568in the variable volume case (\texttt{vvl?} defined), 
     1569where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum}) 
     1570 
     1571$\bullet$ vector invariant form or linear free surface 
     1572(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined): 
     1573\[ 
     1574  % \label{eq:DYN_nxt_vec} 
     1575  \left\{ 
     1576    \begin{aligned} 
     1577      &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\ 
     1578      &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right] 
     1579    \end{aligned} 
     1580  \right. 
     1581\] 
     1582 
     1583$\bullet$ flux form and nonlinear free surface 
     1584(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined): 
     1585\[ 
     1586  % \label{eq:DYN_nxt_flux} 
     1587  \left\{ 
     1588    \begin{aligned} 
     1589      &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\ 
     1590      &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t 
     1591      +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\rdt} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\rdt}} \right] 
     1592    \end{aligned} 
     1593  \right. 
     1594\] 
     1595where RHS is the right hand side of the momentum equation, 
     1596the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient. 
     1597$\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter). 
     1598Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}. 
     1599In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for 
     1600the momentum equations. 
     1601 
     1602Note that with the filtered free surface, 
     1603the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module, 
     1604and only array swapping and Asselin filtering is done in \mdl{dynnxt}. 
     1605 
     1606\onlyinsubfile{\input{../../global/epilogue}} 
     1607 
     1608\end{document} 
Note: See TracChangeset for help on using the changeset viewer.