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.
chap_DYN.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_DYN.tex @ 11582

Last change on this file since 11582 was 11582, checked in by nicolasmartin, 5 years ago

New LaTeX commands \nam and \np to mention namelist content (step 2)
Finally convert \forcode{...} following \np{}{} into optional arg of the new command \np[]{}{}

File size: 92.7 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\onlyinsubfile{\makeindex}
4
5\begin{document}
6% ================================================================
7% Chapter ——— Ocean Dynamics (DYN)
8% ================================================================
9\chapter{Ocean Dynamics (DYN)}
10\label{chap:DYN}
11
12\chaptertoc
13
14Using the representation described in \autoref{chap:DOM},
15several semi-discrete space forms of the dynamical equations are available depending on
16the vertical coordinate used and on the conservation properties of the vorticity term.
17In all the equations presented here, the masking has been omitted for simplicity.
18One must be aware that all the quantities are masked fields and
19that each time an average or difference operator is used, the resulting field is multiplied by a mask.
20
21The prognostic ocean dynamics equation can be summarized as follows:
22\[
23  \text{NXT} = \dbinom  {\text{VOR} + \text{KEG} + \text {ZAD} }
24  {\text{COR} + \text{ADV}                       }
25  + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
26\]
27NXT stands for next, referring to the time-stepping.
28The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that
29are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and
30a vertical advection part (ZAD) in the vector invariant formulation,
31or a Coriolis and advection part (COR+ADV) in the flux formulation.
32The terms following these are the pressure gradient contributions
33(HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient);
34and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF),
35which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules.
36The vertical diffusion term includes the surface and bottom stresses.
37The external forcings and parameterisations require complex inputs
38(surface wind stress calculation using bulk formulae, estimation of mixing coefficients)
39that are carried out in modules SBC, LDF and ZDF and are described in
40\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.
41
42In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence,
43curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module).
44
45The different options available to the user are managed by namelist variables.
46For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx},
47where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
48%If a CPP key is used for this term its name is \key{ttt}.
49The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory,
50and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine.
51
52The user has the option of extracting and outputting each tendency term from the 3D momentum equations
53(\texttt{trddyn?} defined), as described in \autoref{chap:MISC}.
54Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined)
55can be derived from the 3D terms.
56%%%
57\gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does
58MISC correspond to "extracting tendency terms" or "vorticity balance"?}
59
60% ================================================================
61% Sea Surface Height evolution & Diagnostics variables
62% ================================================================
63\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)}
64\label{sec:DYN_divcur_wzv}
65
66%--------------------------------------------------------------------------------------------------------------
67%           Horizontal divergence and relative vorticity
68%--------------------------------------------------------------------------------------------------------------
69\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
70\label{subsec:DYN_divcur}
71
72The vorticity is defined at an $f$-point (\ie\ corner point) as follows:
73\begin{equation}
74  \label{eq:DYN_divcur_cur}
75  \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right]
76      -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right)
77\end{equation}
78
79The horizontal divergence is defined at a $T$-point.
80It is given by:
81\[
82  % \label{eq:DYN_divcur_div}
83  \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
84  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right]
85      +\delta_j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right)
86\]
87
88Note that although the vorticity has the same discrete expression in $z$- and $s$-coordinates,
89its physical meaning is not identical.
90$\zeta$ is a pseudo vorticity along $s$-surfaces
91(only pseudo because $(u,v)$ are still defined along geopotential surfaces,
92but are not necessarily defined at the same depth).
93
94The vorticity and divergence at the \textit{before} step are used in the computation of
95the horizontal diffusion of momentum.
96Note that because they have been calculated prior to the Asselin filtering of the \textit{before} velocities,
97the \textit{before} vorticity and divergence arrays must be included in the restart file to
98ensure perfect restartability.
99The vorticity and divergence at the \textit{now} time step are used for the computation of
100the nonlinear advection and of the vertical velocity respectively.
101
102%--------------------------------------------------------------------------------------------------------------
103%           Sea Surface Height evolution
104%--------------------------------------------------------------------------------------------------------------
105\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
106\label{subsec:DYN_sshwzv}
107
108The sea surface height is given by:
109\begin{equation}
110  \label{eq:DYN_spg_ssh}
111  \begin{aligned}
112    \frac{\partial \eta }{\partial t}
113    &\equiv    \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{  \delta_i \left[ {e_{2u}\,e_{3u}\;u} \right]
114        +\delta_j \left[ {e_{1v}\,e_{3v}\;v} \right]  \right\} }
115    -    \frac{\textit{emp}}{\rho_w }   \\
116    &\equiv    \sum\limits_k {\chi \ e_{3t}}  -  \frac{\textit{emp}}{\rho_w }
117  \end{aligned}
118\end{equation}
119where \textit{emp} is the surface freshwater budget (evaporation minus precipitation),
120expressed in Kg/m$^2$/s (which is equal to mm/s),
121and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation).
122If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then
123\textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.
124The sea-surface height is evaluated using exactly the same time stepping scheme as
125the tracer equation \autoref{eq:TRA_nxt}:
126a leapfrog scheme in combination with an Asselin time filter,
127\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity).
128This is of paramount importance.
129Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
130the sea surface height equation otherwise tracer content will not be conserved
131\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}.
132
133The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
134taking into account the change of the thickness of the levels:
135\begin{equation}
136  \label{eq:DYN_wzv}
137  \left\{
138    \begin{aligned}
139      &\left. w \right|_{k_b-1/2} \quad= 0    \qquad \text{where } k_b \text{ is the level just above the sea floor }   \\
140      &\left. w \right|_{k+1/2}     = \left. w \right|_{k-1/2}  +  \left. e_{3t} \right|_{k}\;  \left. \chi \right|_k
141      - \frac{1} {2 \rdt} \left\left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right)
142    \end{aligned}
143  \right.
144\end{equation}
145
146In the case of a non-linear free surface (\texttt{vvl?}), the top vertical velocity is $-\textit{emp}/\rho_w$,
147as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses,
148re-orientated downward.
149\gmcomment{not sure of this...  to be modified with the change in emp setting}
150In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears.
151The upper boundary condition applies at a fixed level $z=0$.
152The top vertical velocity is thus equal to the divergence of the barotropic transport
153(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}).
154
155Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
156its physical meaning is not the same:
157in the second case, $w$ is the velocity normal to the $s$-surfaces.
158Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to
159the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv}
160(see \autoref{subsec:DOM_Num_Index_vertical}).
161
162
163% ================================================================
164% Coriolis and Advection terms: vector invariant form
165% ================================================================
166\section{Coriolis and advection: vector invariant form}
167\label{sec:DYN_adv_cor_vect}
168%-----------------------------------------nam_dynadv----------------------------------------------------
169
170\begin{listing}
171  \nlst{namdyn_adv}
172  \caption{\forcode{&namdyn_adv}}
173  \label{lst:namdyn_adv}
174\end{listing}
175%-------------------------------------------------------------------------------------------------------------
176
177The vector invariant form of the momentum equations is the one most often used in
178applications of the \NEMO\ ocean model.
179The flux form option (see next section) has been present since version $2$.
180Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and
181momentum advection terms are evaluated using a leapfrog scheme,
182\ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity).
183At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following
184\autoref{chap:LBC}.
185
186% -------------------------------------------------------------------------------------------------------------
187%        Vorticity term
188% -------------------------------------------------------------------------------------------------------------
189\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
190\label{subsec:DYN_vor}
191%------------------------------------------nam_dynvor----------------------------------------------------
192
193\begin{listing}
194  \nlst{namdyn_vor}
195  \caption{\forcode{&namdyn_vor}}
196  \label{lst:namdyn_vor}
197\end{listing}
198%-------------------------------------------------------------------------------------------------------------
199
200Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables.
201Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available:
202conserving potential enstrophy of horizontally non-divergent flow (ENS scheme);
203conserving horizontal kinetic energy (ENE scheme);
204conserving potential enstrophy for the relative vorticity term and
205horizontal kinetic energy for the planetary vorticity term (MIX scheme);
206or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy
207(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}).
208In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
209vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}).
210The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module.
211
212%-------------------------------------------------------------
213%                 enstrophy conserving scheme
214%-------------------------------------------------------------
215\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})}
216\label{subsec:DYN_vor_ens}
217
218In the enstrophy conserving case (ENS scheme),
219the discrete formulation of the vorticity term provides a global conservation of the enstrophy
220($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie\ $\chi$=$0$),
221but does not conserve the total kinetic energy.
222It is given by:
223\begin{equation}
224  \label{eq:DYN_vor_ens}
225  \left\{
226    \begin{aligned}
227      {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i}
228      & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2}    \\
229      {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}
230      & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j}
231    \end{aligned}
232  \right.
233\end{equation}
234
235%-------------------------------------------------------------
236%                 energy conserving scheme
237%-------------------------------------------------------------
238\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})}
239\label{subsec:DYN_vor_ene}
240
241The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy.
242It is given by:
243\begin{equation}
244  \label{eq:DYN_vor_ene}
245  \left\{
246    \begin{aligned}
247      {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
248            \;  \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
249      {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
250            \;  \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} }
251    \end{aligned}
252  \right.
253\end{equation}
254
255%-------------------------------------------------------------
256%                 mix energy/enstrophy conserving scheme
257%-------------------------------------------------------------
258\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})}
259\label{subsec:DYN_vor_mix}
260
261For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used.
262It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term,
263and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term.
264\[
265  % \label{eq:DYN_vor_mix}
266  \left\{ {
267      \begin{aligned}
268        {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}
269          \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} }
270          \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
271              \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
272        {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j
273          \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} }
274          \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
275              \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill
276      \end{aligned}
277    } \right.
278\]
279
280%-------------------------------------------------------------
281%                 energy and enstrophy conserving scheme
282%-------------------------------------------------------------
283\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})}
284\label{subsec:DYN_vor_een}
285
286In both the ENS and ENE schemes,
287it is apparent that the combination of $i$ and $j$ averages of the velocity allows for
288the presence of grid point oscillation structures that will be invisible to the operator.
289These structures are \textit{computational modes} that will be at least partly damped by
290the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term.
291The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field.
292Such noise would result in more noise in the vertical velocity field, an undesirable feature.
293This is a well-known characteristic of $C$-grid discretization where
294$u$ and $v$ are located at different grid points,
295a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid.
296\gmcomment{ To circumvent this, Adcroft (ADD REF HERE)
297Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
298
299A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}.
300The idea is to get rid of the double averaging by considering triad combinations of vorticity.
301It is noteworthy that this solution is conceptually quite similar to the one proposed by
302\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}).
303
304The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
305for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme.
306First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
307\[
308  % \label{eq:DYN_pot_vor}
309  q  = \frac{\zeta +f} {e_{3f} }
310\]
311where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}),
312the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
313\begin{equation}
314  \label{eq:DYN_een_e3f}
315  e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
316\end{equation}
317
318%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
319\begin{figure}[!ht]
320  \centering
321  \includegraphics[width=0.66\textwidth]{Fig_DYN_een_triad}
322  \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{
323    Triads used in the energy and enstrophy conserving scheme (EEN) for
324    $u$-component (upper panel) and $v$-component (lower panel).}
325  \label{fig:DYN_een_triad}
326\end{figure}
327% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
328
329A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
330It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks
331(\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}).
332The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
333extends by continuity the value of $e_{3f}$ into the land areas.
334This case introduces a sub-grid-scale topography at f-points
335(with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry)
336that tends to reinforce the topostrophy of the flow
337(\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}.
338
339Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
340the following triad combinations of the neighbouring potential vorticities defined at f-points
341(\autoref{fig:DYN_een_triad}):
342\begin{equation}
343  \label{eq:DYN_Q_triads}
344  _i^j \mathbb{Q}^{i_p}_{j_p}
345  = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right)
346\end{equation}
347where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$.
348
349Finally, the vorticity terms are represented as:
350\begin{equation}
351  \label{eq:DYN_vor_een}
352  \left\{ {
353      \begin{aligned}
354        +q\,e_3 \, v    &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}}
355        {^{i+1/2-i_p}_j}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{1v}\,e_{3v} \;\right)^{i+1/2-i_p}_{j+j_p}   \\
356        - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}}
357        {^i_{j+1/2-j_p}}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{2u}\,e_{3u} \;\right)^{i+i_p}_{j+1/2-j_p}   \\
358      \end{aligned}
359    } \right.
360\end{equation}
361
362This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
363It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow
364(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}).
365Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
366the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
367Furthermore, used in combination with a partial steps representation of bottom topography,
368it improves the interaction between current and topography,
369leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.
370
371%--------------------------------------------------------------------------------------------------------------
372%           Kinetic Energy Gradient term
373%--------------------------------------------------------------------------------------------------------------
374\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})}
375\label{subsec:DYN_keg}
376
377As demonstrated in \autoref{apdx:INVARIANTS},
378there is a single discrete formulation of the kinetic energy gradient term that,
379together with the formulation chosen for the vertical advection (see below),
380conserves the total kinetic energy:
381\[
382  % \label{eq:DYN_keg}
383  \left\{
384    \begin{aligned}
385      -\frac{1}{2 \; e_{1u} }  & \ \delta_{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\
386      -\frac{1}{2 \; e_{2v} }  & \ \delta_{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]
387    \end{aligned}
388  \right.
389\]
390
391%--------------------------------------------------------------------------------------------------------------
392%           Vertical advection term
393%--------------------------------------------------------------------------------------------------------------
394\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})}
395\label{subsec:DYN_zad}
396
397The discrete formulation of the vertical advection, t
398ogether with the formulation chosen for the gradient of kinetic energy (KE) term,
399conserves the total kinetic energy.
400Indeed, the change of KE due to the vertical advection is exactly balanced by
401the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}).
402\[
403  % \label{eq:DYN_zad}
404  \left\{
405    \begin{aligned}
406      -\frac{1} {e_{1u}\,e_{2u}\,e_{3u}} &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,i+1/2\;\delta_{k+1/2} \left[ u \right]\  }^{\,k}  \\
407      -\frac{1} {e_{1v}\,e_{2v}\,e_{3v}}  &\ \overline{\ \overline{ e_{1t}\,e_{2t}\;w } ^{\,j+1/2\;\delta_{k+1/2} \left[ u \right]\  }^{\,k}
408    \end{aligned}
409  \right.
410\]
411When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts},
412a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term.
413This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
414Note that in this case,
415a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability,
416an 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}).
417
418
419% ================================================================
420% Coriolis and Advection : flux form
421% ================================================================
422\section{Coriolis and advection: flux form}
423\label{sec:DYN_adv_cor_flux}
424%------------------------------------------nam_dynadv----------------------------------------------------
425
426%-------------------------------------------------------------------------------------------------------------
427
428Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables.
429In the flux form (as in the vector invariant form),
430the Coriolis and momentum advection terms are evaluated using a leapfrog scheme,
431\ie\ the velocity appearing in their expressions is centred in time (\textit{now} velocity).
432At the lateral boundaries either free slip,
433no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}.
434
435
436%--------------------------------------------------------------------------------------------------------------
437%           Coriolis plus curvature metric terms
438%--------------------------------------------------------------------------------------------------------------
439\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
440\label{subsec:DYN_cor_flux}
441
442In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term.
443This altered Coriolis parameter is thus discretised at $f$-points.
444It is given by:
445\begin{multline*}
446  % \label{eq:DYN_cor_metric}
447  f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right\\
448  \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right]
449      -  \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right]  }  \ \right)
450\end{multline*}
451
452Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to
453compute the product of the Coriolis parameter and the vorticity.
454However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date.
455This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity).
456
457%--------------------------------------------------------------------------------------------------------------
458%           Flux form Advection term
459%--------------------------------------------------------------------------------------------------------------
460\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})}
461\label{subsec:DYN_adv_flux}
462
463The discrete expression of the advection term is given by:
464\[
465  % \label{eq:DYN_adv}
466  \left\{
467    \begin{aligned}
468      \frac{1}{e_{1u}\,e_{2u}\,e_{3u}}
469      \left(      \delta_{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i       }  \ u_t      \right]
470        + \delta_{j       } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2\ u_f      \right] \right\ \;   \\
471      \left.   + \delta_{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2\ u_{uw} \right] \right)   \\
472      \\
473      \frac{1}{e_{1v}\,e_{2v}\,e_{3v}}
474      \left(     \delta_{i       } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f       \right]
475        + \delta_{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i       } \ v_t       \right] \right\ \, \, \\
476      \left+ \delta_{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw}  \right] \right) \\
477    \end{aligned}
478  \right.
479\]
480
481Two advection schemes are available:
482a $2^{nd}$ order centered finite difference scheme, CEN2,
483or a $3^{rd}$ order upstream biased scheme, UBS.
484The latter is described in \citet{shchepetkin.mcwilliams_OM05}.
485The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_ubs}{ln\_dynadv\_ubs}.
486In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
487$u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-,
488and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$.
489
490%-------------------------------------------------------------
491%                 2nd order centred scheme
492%-------------------------------------------------------------
493\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})}
494\label{subsec:DYN_adv_cen2}
495
496In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points:
497\begin{equation}
498  \label{eq:DYN_adv_cen2}
499  \left\{
500    \begin{aligned}
501      u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\
502      v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j    \quad &  v_{vw}^{cen2} &=\overline v ^{k+1/2}  \\
503    \end{aligned}
504  \right.
505\end{equation}
506
507The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema).
508It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to
509produce a sensible solution.
510The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter,
511so $u$ and $v$ are the \emph{now} velocities.
512
513%-------------------------------------------------------------
514%                 UBS scheme
515%-------------------------------------------------------------
516\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})}
517\label{subsec:DYN_adv_ubs}
518
519The UBS advection scheme is an upstream biased third order scheme based on
520an upstream-biased parabolic interpolation.
521For example, the evaluation of $u_T^{ubs} $ is done as follows:
522\begin{equation}
523  \label{eq:DYN_adv_ubs}
524  u_T^{ubs} =\overline u ^i-\;\frac{1}{6}
525  \begin{cases}
526    u"_{i-1/2}&   \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  \geqslant 0$ }    \\
527    u"_{i+1/2}&   \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  < 0$ }
528  \end{cases}
529\end{equation}
530where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$.
531This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error
532\citep{shchepetkin.mcwilliams_OM05}.
533The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}.
534It is a relatively good compromise between accuracy and smoothness.
535It is not a \emph{positive} scheme, meaning that false extrema are permitted.
536But the amplitudes of the false extrema are significantly reduced over those in the centred second order method.
537As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum
538(\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}),
539and it is recommended to do so.
540
541The UBS scheme is not used in all directions.
542In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and
543$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used.
544UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the
545sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }
546
547For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}),
548which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time),
549while the second term, which is the diffusion part of the scheme,
550is evaluated using the \textit{before} velocity (forward in time).
551This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
552
553Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
554one coefficient.
555Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
556This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
557Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
558
559Note also that in the current version of \mdl{dynadv\_ubs},
560there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS.
561This is an error and should be suppressed soon.
562%%%
563\gmcomment{action :  this have to be done}
564%%%
565
566% ================================================================
567%           Hydrostatic pressure gradient term
568% ================================================================
569\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
570\label{sec:DYN_hpg}
571%------------------------------------------nam_dynhpg---------------------------------------------------
572
573\begin{listing}
574  \nlst{namdyn_hpg}
575  \caption{\forcode{&namdyn_hpg}}
576  \label{lst:namdyn_hpg}
577\end{listing}
578%-------------------------------------------------------------------------------------------------------------
579
580Options are defined through the \nam{dyn_hpg}{dyn\_hpg} namelist variables.
581The key distinction between the different algorithms used for
582the hydrostatic pressure gradient is the vertical coordinate used,
583since HPG is a \emph{horizontal} pressure gradient, \ie\ computed along geopotential surfaces.
584As a result, any tilt of the surface of the computational levels will require a specific treatment to
585compute the hydrostatic pressure gradient.
586
587The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme,
588\ie\ the density appearing in its expression is centred in time (\emph{now} $\rho$),
589or a semi-implcit scheme.
590At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied.
591
592%--------------------------------------------------------------------------------------------------------------
593%           z-coordinate with full step
594%--------------------------------------------------------------------------------------------------------------
595\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco})]{Full step $Z$-coordinate (\protect\np{ln_dynhpg_zco}{ln\_dynhpg\_zco})}
596\label{subsec:DYN_hpg_zco}
597
598The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface.
599However, the pressure is large at great depth while its horizontal gradient is several orders of magnitude smaller.
600This may lead to large truncation errors in the pressure gradient terms.
601Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows:
602
603for $k=km$ (surface layer, $jk=1$ in the code)
604\begin{equation}
605  \label{eq:DYN_hpg_zco_surf}
606  \left\{
607    \begin{aligned}
608      \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k=km}
609      &= \frac{1}{2} g \   \left. \delta_{i+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
610      \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k=km}
611      &= \frac{1}{2} g \   \left. \delta_{j+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
612    \end{aligned}
613  \right.
614\end{equation}
615
616for $1<k<km$ (interior layer)
617\begin{equation}
618  \label{eq:DYN_hpg_zco}
619  \left\{
620    \begin{aligned}
621      \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k}
622      &=             \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k-1}
623      +    \frac{1}{2}\;g\;   \left. \delta_{i+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
624      \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k}
625      &=                \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k-1}
626      +    \frac{1}{2}\;g\;   \left. \delta_{j+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
627    \end{aligned}
628  \right.
629\end{equation}
630
631Note that the $1/2$ factor in (\autoref{eq:DYN_hpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
632the vertical derivative of the scale factor at the surface level ($z=0$).
633Note also that in case of variable volume level (\texttt{vvl?} defined),
634the surface pressure gradient is included in \autoref{eq:DYN_hpg_zco_surf} and
635\autoref{eq:DYN_hpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
636
637%--------------------------------------------------------------------------------------------------------------
638%           z-coordinate with partial step
639%--------------------------------------------------------------------------------------------------------------
640\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps})]{Partial step $Z$-coordinate (\protect\np{ln_dynhpg_zps}{ln\_dynhpg\_zps})}
641\label{subsec:DYN_hpg_zps}
642
643With partial bottom cells, tracers in horizontally adjacent cells generally live at different depths.
644Before taking horizontal gradients between these tracer points,
645a linear interpolation is used to approximate the deeper tracer as if
646it actually lived at the depth of the shallower tracer point.
647
648Apart from this modification,
649the horizontal hydrostatic pressure gradient evaluated in the $z$-coordinate with partial step is exactly as in
650the pure $z$-coordinate case.
651As explained in detail in section \autoref{sec:TRA_zpshde},
652the nonlinearity of pressure effects in the equation of state is such that
653it is better to interpolate temperature and salinity vertically before computing the density.
654Horizontal gradients of temperature and salinity are needed for the TRA modules,
655which is the reason why the horizontal gradients of density at the deepest model level are computed in
656module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}.
657
658%--------------------------------------------------------------------------------------------------------------
659%           s- and s-z-coordinates
660%--------------------------------------------------------------------------------------------------------------
661\subsection{$S$- and $Z$-$S$-coordinates}
662\label{subsec:DYN_hpg_sco}
663
664Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers
665(\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}).
666A number of different pressure gradient options are coded but the ROMS-like,
667density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
668
669$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco})
670\begin{equation}
671  \label{eq:DYN_hpg_sco}
672  \left\{
673    \begin{aligned}
674      - \frac{1}                 {\rho_o \, e_{1u}} \;   \delta_{i+1/2} \left[  p^h  \right]
675      + \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta_{i+1/2} \left[  z_t   \right]    \\
676      - \frac{1}                 {\rho_o \, e_{2v}} \;   \delta_{j+1/2} \left[  p^h  \right]
677      + \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta_{j+1/2} \left[  z_t   \right]    \\
678    \end{aligned}
679  \right.
680\end{equation}
681
682Where the first term is the pressure gradient along coordinates,
683computed as in \autoref{eq:DYN_hpg_zco_surf} - \autoref{eq:DYN_hpg_zco},
684and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
685($e_{3w}$).
686
687$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).
688This scheme need the activation of ice shelf cavities (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
689
690$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj})
691
692$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}
693(\np[=.true.]{ln_dynhpg_djc}{ln\_dynhpg\_djc}) (currently disabled; under development)
694
695Note that expression \autoref{eq:DYN_hpg_sco} is commonly used when the variable volume formulation is activated
696(\texttt{vvl?}) because in that case, even with a flat bottom,
697the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
698The pressure jacobian scheme (\np[=.true.]{ln_dynhpg_prj}{ln\_dynhpg\_prj}) is available as
699an improved option to \np[=.true.]{ln_dynhpg_sco}{ln\_dynhpg\_sco} when \texttt{vvl?} is active.
700The pressure Jacobian scheme uses a constrained cubic spline to
701reconstruct the density profile across the water column.
702This method maintains the monotonicity between the density nodes.
703The pressure can be calculated by analytical integration of the density profile and
704a pressure Jacobian method is used to solve the horizontal pressure gradient.
705This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.
706
707\subsection{Ice shelf cavity}
708\label{subsec:DYN_hpg_isf}
709
710Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and
711the pressure gradient due to the ocean load (\np[=.true.]{ln_dynhpg_isf}{ln\_dynhpg\_isf}).\\
712
713The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium.
714The top pressure is computed integrating from surface to the base of the ice shelf a reference density profile
715(prescribed as density of a water at 34.4 PSU and -1.9\deg{C}) and
716corresponds to the water replaced by the ice shelf.
717This top pressure is constant over time.
718A detailed description of this method is described in \citet{losch_JGR08}.\\
719
720The pressure gradient due to ocean load is computed using the expression \autoref{eq:DYN_hpg_sco} described in
721\autoref{subsec:DYN_hpg_sco}.
722
723%--------------------------------------------------------------------------------------------------------------
724%           Time-scheme
725%--------------------------------------------------------------------------------------------------------------
726\subsection[Time-scheme (\forcode{ln_dynhpg_imp})]{Time-scheme (\protect\np{ln_dynhpg_imp}{ln\_dynhpg\_imp})}
727\label{subsec:DYN_hpg_imp}
728
729The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and
730therefore the density used in all discrete expressions given above is the  \textit{now} density,
731computed from the \textit{now} temperature and salinity.
732In some specific cases
733(usually high resolution simulations over an ocean domain which includes weakly stratified regions)
734the physical phenomenon that controls the time-step is internal gravity waves (IGWs).
735A semi-implicit scheme for doubling the stability limit associated with IGWs can be used
736\citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}.
737It involves the evaluation of the hydrostatic pressure gradient as
738an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$
739(\ie\ \textit{before}, \textit{now} and  \textit{after} time-steps),
740rather than at the central time level $t$ only, as in the standard leapfrog scheme.
741
742$\bullet$ leapfrog scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
743
744\begin{equation}
745  \label{eq:DYN_hpg_lf}
746  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
747  -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right]
748\end{equation}
749
750$\bullet$ semi-implicit scheme (\np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}):
751\begin{equation}
752  \label{eq:DYN_hpg_imp}
753  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
754  -\frac{1}{4\,\rho_o \,e_{1u} } \delta_{i+1/2} \left[ p_h^{t+\rdt} +2\,p_h^t +p_h^{t-\rdt}  \right]
755\end{equation}
756
757The semi-implicit time scheme \autoref{eq:DYN_hpg_imp} is made possible without
758significant additional computation since the density can be updated to time level $t+\rdt$ before
759computing the horizontal hydrostatic pressure gradient.
760It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using
761\autoref{eq:DYN_hpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:DYN_hpg_lf}.
762Note that \autoref{eq:DYN_hpg_imp} is equivalent to applying a time filter to the pressure gradient to
763eliminate high frequency IGWs.
764Obviously, when using \autoref{eq:DYN_hpg_imp},
765the doubling of the time-step is achievable only if no other factors control the time-step,
766such as the stability limits associated with advection or diffusion.
767
768In practice, the semi-implicit scheme is used when \np[=.true.]{ln_dynhpg_imp}{ln\_dynhpg\_imp}.
769In this case, we choose to apply the time filter to temperature and salinity used in the equation of state,
770instead of applying it to the hydrostatic pressure or to the density,
771so that no additional storage array has to be defined.
772The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:
773\[
774  % \label{eq:DYN_rho_flt}
775  \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
776  \quad    \text{with}  \quad
777  \widetilde{X} = 1 / 4 \left(  X^{t+\rdt} +2 \,X^t + X^{t-\rdt\right)
778\]
779
780Note that in the semi-implicit case, it is necessary to save the filtered density,
781an extra three-dimensional field, in the restart file to restart the model with exact reproducibility.
782This option is controlled by  \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter.
783
784% ================================================================
785% Surface Pressure Gradient
786% ================================================================
787\section[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})}
788\label{sec:DYN_spg}
789%-----------------------------------------nam_dynspg----------------------------------------------------
790
791\begin{listing}
792  \nlst{namdyn_spg}
793  \caption{\forcode{&namdyn_spg}}
794  \label{lst:namdyn_spg}
795\end{listing}
796%------------------------------------------------------------------------------------------------------------
797
798Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.
799The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:MB_hor_pg}).
800The main distinction is between the fixed volume case (linear free surface) and
801the variable volume case (nonlinear free surface, \texttt{vvl?} is defined).
802In the linear free surface case (\autoref{subsec:MB_free_surface})
803the vertical scale factors $e_{3}$ are fixed in time,
804while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
805With both linear and nonlinear free surface, external gravity waves are allowed in the equations,
806which imposes a very small time step when an explicit time stepping is used.
807Two methods are proposed to allow a longer time step for the three-dimensional equations:
808the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:MB_flt?}),
809and the split-explicit free surface described below.
810The extra term introduced in the filtered method is calculated implicitly,
811so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
812
813
814The form of the surface pressure gradient term depends on how the user wants to
815handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:MB_hor_pg}).
816Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx):
817an explicit formulation which requires a small time step;
818a filtered free surface formulation which allows a larger time step by
819adding a filtering term into the momentum equation;
820and a split-explicit free surface formulation, described below, which also allows a larger time step.
821
822The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it.
823As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
824
825
826%--------------------------------------------------------------------------------------------------------------
827% Explicit free surface formulation
828%--------------------------------------------------------------------------------------------------------------
829\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
830\label{subsec:DYN_spg_exp}
831
832In the explicit free surface formulation (\np{ln_dynspg_exp}{ln\_dynspg\_exp} set to true),
833the model time step is chosen to be small enough to resolve the external gravity waves
834(typically a few tens of seconds).
835The surface pressure gradient, evaluated using a leap-frog scheme (\ie\ centered in time),
836is thus simply given by :
837\begin{equation}
838  \label{eq:DYN_spg_exp}
839  \left\{
840    \begin{aligned}
841      - \frac{1}{e_{1u}\,\rho_o} \; \delta_{i+1/2} \left[  \,\rho \,\eta\,  \right]    \\
842      - \frac{1}{e_{2v}\,\rho_o} \; \delta_{j+1/2} \left[  \,\rho \,\eta\,  \right]
843    \end{aligned}
844  \right.
845\end{equation}
846
847Note that in the non-linear free surface case (\ie\ \texttt{vvl?} defined),
848the surface pressure gradient is already included in the momentum tendency through
849the level thickness variation allowed in the computation of the hydrostatic pressure gradient.
850Thus, nothing is done in the \mdl{dynspg\_exp} module.
851
852%--------------------------------------------------------------------------------------------------------------
853% Split-explict free surface formulation
854%--------------------------------------------------------------------------------------------------------------
855\subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})}
856\label{subsec:DYN_spg_ts}
857%------------------------------------------namsplit-----------------------------------------------------------
858%
859%\nlst{namsplit}
860%-------------------------------------------------------------------------------------------------------------
861
862The split-explicit free surface formulation used in \NEMO\ (\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true),
863also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}.
864The general idea is to solve the free surface equation and the associated barotropic velocity equations with
865a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
866(\autoref{fig:DYN_spg_ts}).
867The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through
868the \np{nn_baro}{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
869This parameter can be optionally defined automatically (\np[=.true.]{ln_bt_nn_auto}{ln\_bt\_nn\_auto}) considering that
870the stability of the barotropic system is essentially controled by external waves propagation.
871Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry.
872Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn_bt_cmax}{rn\_bt\_cmax}.
873
874%%%
875The barotropic mode solves the following equations:
876% \begin{subequations}
877%  \label{eq:DYN_BT}
878\begin{equation}
879  \label{eq:DYN_BT_dyn}
880  \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}=
881  -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h}
882  -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}}
883\end{equation}
884\[
885  % \label{eq:DYN_BT_ssh}
886  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E
887\]
888% \end{subequations}
889where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes,
890surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency.
891The third term on the right hand side of \autoref{eq:DYN_BT_dyn} represents the bottom stress
892(see section \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.
893Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
894detailed in \citet{shchepetkin.mcwilliams_OM05}.
895AB3-AM4 coefficients used in \NEMO\ follow the second-order accurate,
896"multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09}
897(see their figure 12, lower left).
898
899%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
900\begin{figure}[!t]
901  \centering
902  \includegraphics[width=0.66\textwidth]{Fig_DYN_dynspg_ts}
903  \caption[Split-explicit time stepping scheme for the external and internal modes]{
904    Schematic of the split-explicit time stepping scheme for the external and internal modes.
905    Time increases to the right.
906    In this particular exemple,
907    a boxcar averaging window over \np{nn_baro}{nn\_baro} barotropic time steps is used
908    (\np[=1]{nn_bt_flt}{nn\_bt\_flt}) and \np[=5]{nn_baro}{nn\_baro}.
909    Internal mode time steps (which are also the model time steps) are denoted by
910    $t-\rdt$, $t$ and $t+\rdt$.
911    Variables with $k$ superscript refer to instantaneous barotropic variables,
912    $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary
913    (red vertical bars) and secondary weights (blue vertical bars).
914    The former are used to obtain time filtered quantities at $t+\rdt$ while
915    the latter are used to obtain time averaged transports to advect tracers.
916    a) Forward time integration:
917    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}\protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
918    b) Centred time integration:
919    \protect\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}, \protect\np[=.true.]{ln_bt_av}{ln\_bt\_av}.
920    c) Forward time integration with no time filtering (POM-like scheme):
921    \protect\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}\protect\np[=.false.]{ln_bt_av}{ln\_bt\_av}.}
922  \label{fig:DYN_spg_ts}
923\end{figure}
924%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
925
926In the default case (\np[=.true.]{ln_bt_fw}{ln\_bt\_fw}),
927the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps
928(\autoref{fig:DYN_spg_ts}a).
929To avoid aliasing of fast barotropic motions into three dimensional equations,
930time filtering is eventually applied on barotropic quantities (\np[=.true.]{ln_bt_av}{ln\_bt\_av}).
931In that case, the integration is extended slightly beyond \textit{after} time step to
932provide time filtered quantities.
933These are used for the subsequent initialization of the barotropic mode in the following baroclinic step.
934Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme,
935asselin filtering is not applied to barotropic quantities.\\
936Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step
937(\np[=.false.]{ln_bt_fw}{ln\_bt\_fw}).
938Although more computationaly expensive ( \np{nn_baro}{nn\_baro} additional iterations are indeed necessary),
939the baroclinic to barotropic forcing term given at \textit{now} time step become centred in
940the middle of the integration window.
941It can easily be shown that this property removes part of splitting errors between modes,
942which increases the overall numerical robustness.
943%references to Patrick Marsaleix' work here. Also work done by SHOM group.
944
945%%%
946
947As far as tracer conservation is concerned,
948barotropic velocities used to advect tracers must also be updated at \textit{now} time step.
949This implies to change the traditional order of computations in \NEMO:
950most of momentum trends (including the barotropic mode calculation) updated first, tracers' after.
951This \textit{de facto} makes semi-implicit hydrostatic pressure gradient
952(see section \autoref{subsec:DYN_hpg_imp})
953and time splitting not compatible.
954Advective barotropic velocities are obtained by using a secondary set of filtering weights,
955uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}).
956Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to
957obtain exact conservation.
958
959%%%
960
961One can eventually choose to feedback instantaneous values by not using any time filter
962(\np[=.false.]{ln_bt_av}{ln\_bt\_av}).
963In that case, external mode equations are continuous in time,
964\ie\ they are not re-initialized when starting a new sub-stepping sequence.
965This is the method used so far in the POM model, the stability being maintained by
966refreshing at (almost) each barotropic time step advection and horizontal diffusion terms.
967Since the latter terms have not been added in \NEMO\ for computational efficiency,
968removing time filtering is not recommended except for debugging purposes.
969This may be used for instance to appreciate the damping effect of the standard formulation on
970external gravity waves in idealized or weakly non-linear cases.
971Although the damping is lower than for the filtered free surface,
972it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave.
973
974%>>>>>===============
975\gmcomment{               %%% copy from griffies Book
976
977\textbf{title: Time stepping the barotropic system }
978
979Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$.
980Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using
981the small barotropic time step $\rdt$.
982We have
983
984\[
985  % \label{eq:DYN_spg_ts_eta}
986  \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
987  = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right]
988\]
989\begin{multline*}
990  % \label{eq:DYN_spg_ts_u}
991  \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}\\
992  = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})
993    - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right]
994\end{multline*}
995\
996
997In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with
998the barotropic time steps.
999The $\tau$ time label on $\eta^{(b)}$ and $U^{(b)}$ denotes the baroclinic time at which
1000the vertically integrated forcing $\textbf{M}(\tau)$
1001(note that this forcing includes the surface freshwater forcing),
1002the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$,
1003and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over
1004a single cycle.
1005This is also the time that sets the barotropic time steps via
1006\[
1007  % \label{eq:DYN_spg_ts_t}
1008  t_n=\tau+n\rdt
1009\]
1010with $n$ an integer.
1011The density scaled surface pressure is evaluated via
1012\[
1013  % \label{eq:DYN_spg_ts_ps}
1014  p_s^{(b)}(\tau,t_{n}) =
1015  \begin{cases}
1016    g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o  &      \text{non-linear case} \\
1017    g \;\eta_s^{(b)}(\tau,t_{n})  &      \text{linear case}
1018  \end{cases}
1019\]
1020To get started, we assume the following initial conditions
1021\[
1022  % \label{eq:DYN_spg_ts_eta}
1023  \begin{split}
1024    \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)}    \\
1025    \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
1026  \end{split}
1027\]
1028with
1029\[
1030  % \label{eq:DYN_spg_ts_etaF}
1031  \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
1032\]
1033the time averaged surface height taken from the previous barotropic cycle.
1034Likewise,
1035\[
1036  % \label{eq:DYN_spg_ts_u}
1037  \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ \\
1038  \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
1039\]
1040with
1041\[
1042  % \label{eq:DYN_spg_ts_u}
1043  \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n})
1044\]
1045the time averaged vertically integrated transport.
1046Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration.
1047
1048Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ ,
1049the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at
1050baroclinic time $\tau + \rdt \tau$
1051\[
1052  % \label{eq:DYN_spg_ts_u}
1053  \textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n})
1054\]
1055The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using
1056the following form
1057
1058\begin{equation}
1059  \label{eq:DYN_spg_ts_ssh}
1060  \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]
1061\end{equation}
1062
1063The use of this "big-leap-frog" scheme for the surface height ensures compatibility between
1064the mass/volume budgets and the tracer budgets.
1065More discussion of this point is provided in Chapter 10 (see in particular Section 10.2).
1066
1067In general, some form of time filter is needed to maintain integrity of the surface height field due to
1068the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
1069We have tried various forms of such filtering,
1070with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
1071its stability and reasonably good maintenance of tracer conservation properties (see ??).
1072
1073\begin{equation}
1074  \label{eq:DYN_spg_ts_sshf}
1075  \eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)}
1076\end{equation}
1077Another approach tried was
1078
1079\[
1080  % \label{eq:DYN_spg_ts_sshf2}
1081  \eta^{F}(\tau-\Delta) = \eta(\tau)
1082  + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
1083    + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
1084\]
1085
1086which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$.
1087This isolation allows for an easy check that tracer conservation is exact when
1088eliminating tracer and surface height time filtering (see ?? for more complete discussion).
1089However, in the general case with a non-zero $\alpha$,
1090the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.
1091
1092}            %%end gm comment (copy of griffies book)
1093
1094%>>>>>===============
1095
1096
1097%--------------------------------------------------------------------------------------------------------------
1098% Filtered free surface formulation
1099%--------------------------------------------------------------------------------------------------------------
1100\subsection{Filtered free surface (\forcode{dynspg_flt?})}
1101\label{subsec:DYN_spg_fltp}
1102
1103The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.
1104The extra term introduced in the equations (see \autoref{subsec:MB_free_surface}) is solved implicitly.
1105The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
1106
1107%% gm %%======>>>>   given here the discrete eqs provided to the solver
1108\gmcomment{               %%% copy from chap-model basics
1109  \[
1110    % \label{eq:DYN_spg_flt}
1111    \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}}
1112    - g \nabla \left( \tilde{\rho} \ \eta \right)
1113    - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)
1114  \]
1115  where $T_c$, is a parameter with dimensions of time which characterizes the force,
1116  $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density,
1117  and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
1118  non-linear and viscous terms in \autoref{eq:MB_dyn}.
1119}   %end gmcomment
1120
1121Note that in the linear free surface formulation (\texttt{vvl?} not defined),
1122the ocean depth is time-independent and so is the matrix to be inverted.
1123It is computed once and for all and applies to all ocean time steps.
1124
1125% ================================================================
1126% Lateral diffusion term
1127% ================================================================
1128\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
1129\label{sec:DYN_ldf}
1130%------------------------------------------nam_dynldf----------------------------------------------------
1131
1132\begin{listing}
1133  \nlst{namdyn_ldf}
1134  \caption{\forcode{&namdyn_ldf}}
1135  \label{lst:namdyn_ldf}
1136\end{listing}
1137%-------------------------------------------------------------------------------------------------------------
1138
1139Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables.
1140The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators.
1141The coefficients may be constant or spatially variable;
1142the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}).
1143The lateral diffusion of momentum is evaluated using a forward scheme,
1144\ie\ the velocity appearing in its expression is the \textit{before} velocity in time,
1145except for the pure vertical component that appears when a tensor of rotation is used.
1146This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}).
1147
1148At the lateral boundaries either free slip,
1149no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}).
1150
1151\gmcomment{
1152  Hyperviscous operators are frequently used in the simulation of turbulent flows to
1153  control the dissipation of unresolved small scale features.
1154  Their primary role is to provide strong dissipation at the smallest scale supported by
1155  the grid while minimizing the impact on the larger scale features.
1156  Hyperviscous operators are thus designed to be more scale selective than the traditional,
1157  physically motivated Laplace operator.
1158  In finite difference methods,
1159  the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since
1160  its damping time (\ie\ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$
1161  (so that short waves damped more rapidelly than long ones),
1162  whereas the Laplace operator damping time scales only like $\lambda^{-2}$.
1163}
1164
1165% ================================================================
1166\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
1167\label{subsec:DYN_ldf_lap}
1168
1169For lateral iso-level diffusion, the discrete operator is:
1170\begin{equation}
1171  \label{eq:DYN_ldf_lap}
1172  \left\{
1173    \begin{aligned}
1174      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
1175          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
1176        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
1177      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
1178          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
1179        {A_f^{lm} \;e_{3f} \zeta } \right]
1180    \end{aligned}
1181  \right.
1182\end{equation}
1183
1184As explained in \autoref{subsec:MB_ldf},
1185this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and
1186ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.
1187
1188%--------------------------------------------------------------------------------------------------------------
1189%           Rotated laplacian operator
1190%--------------------------------------------------------------------------------------------------------------
1191\subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})}
1192\label{subsec:DYN_ldf_iso}
1193
1194A rotation of the lateral momentum diffusion operator is needed in several cases:
1195for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and
1196for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or
1197geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate.
1198In the partial step case, coordinates are horizontal except at the deepest level and
1199no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
1200The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
1201each momentum component.
1202It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.
1203The resulting discrete representation is:
1204\begin{equation}
1205  \label{eq:DYN_ldf_iso}
1206  \begin{split}
1207    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
1208    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_T^{lm}  \left(
1209              {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]
1210                -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}
1211            \right)} \right]}    \right. \\
1212    & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}
1213            }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}
1214            \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}
1215        \right)} \right] \\
1216    &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline
1217              {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }
1218        \right.} \right. \\
1219    &  \ \qquad \qquad \qquad \quad\
1220    - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\
1221    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1222                +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}
1223                \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\
1224    D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\
1225    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_f^{lm}  \left(
1226              {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]
1227                -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}
1228            \right)} \right]}    \right. \\
1229    & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}
1230            }\,\delta_{j} [v] - e_{1t}\, r_{2t}
1231            \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}
1232        \right)} \right] \\
1233    & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline
1234              {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\
1235    &  \ \qquad \qquad \qquad \quad\
1236    - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\
1237    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1238                +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}
1239                \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\}
1240  \end{split}
1241\end{equation}
1242where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and
1243the surface of computation ($z$- or $s$-surfaces).
1244The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).
1245
1246%--------------------------------------------------------------------------------------------------------------
1247%           Iso-level bilaplacian operator
1248%--------------------------------------------------------------------------------------------------------------
1249\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})}
1250\label{subsec:DYN_ldf_bilap}
1251
1252The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.
1253It requires an additional assumption on boundary conditions:
1254the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
1255while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).
1256%%%
1257\gmcomment{add a remark on the the change in the position of the coefficient}
1258%%%
1259
1260% ================================================================
1261%           Vertical diffusion term
1262% ================================================================
1263\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
1264\label{sec:DYN_zdf}
1265%----------------------------------------------namzdf------------------------------------------------------
1266
1267%-------------------------------------------------------------------------------------------------------------
1268
1269Options are defined through the \nam{zdf}{zdf} namelist variables.
1270The 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.
1271Two time stepping schemes can be used for the vertical diffusion term:
1272$(a)$ a forward time differencing scheme
1273(\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or
1274$(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp})
1275(see \autoref{chap:TD}).
1276Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics.
1277
1278The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
1279The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:
1280\[
1281  % \label{eq:DYN_zdf}
1282  \left\{
1283    \begin{aligned}
1284      D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
1285        \ \delta_{k+1/2} [\,u\,]         \right]     \\
1286      \\
1287      D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
1288        \ \delta_{k+1/2} [\,v\,]         \right]
1289    \end{aligned}
1290  \right.
1291\]
1292where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.
1293The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).
1294
1295The surface boundary condition on momentum is the stress exerted by the wind.
1296At the surface, the momentum fluxes are prescribed as the boundary condition on
1297the vertical turbulent momentum fluxes,
1298\begin{equation}
1299  \label{eq:DYN_zdf_sbc}
1300  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
1301  = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
1302\end{equation}
1303where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in
1304the (\textbf{i},\textbf{j}) coordinate system.
1305The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in
1306the vertical over the mixed layer depth.
1307If the vertical mixing coefficient is small (when no mixed layer scheme is used)
1308the surface stress enters only the top model level, as a body force.
1309The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).
1310
1311The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
1312(see \autoref{sec:ZDF_drg})
1313
1314% ================================================================
1315% External Forcing
1316% ================================================================
1317\section{External forcings}
1318\label{sec:DYN_forcing}
1319
1320Besides the surface and bottom stresses (see the above section)
1321which are introduced as boundary conditions on the vertical mixing,
1322three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
1323
1324(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
1325the atmospheric pressure is taken into account when computing the surface pressure gradient.
1326
1327(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_tide}),
1328the tidal potential is taken into account when computing the surface pressure gradient.
1329
1330(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and LIM or CICE is used
1331(\ie\ when the sea-ice is embedded in the ocean),
1332the snow-ice mass is taken into account when computing the surface pressure gradient.
1333
1334
1335\gmcomment{ missing : the lateral boundary condition !!!   another external forcing
1336 }
1337
1338% ================================================================
1339% Wetting and drying
1340% ================================================================
1341\section{Wetting and drying }
1342\label{sec:DYN_wetdry}
1343
1344There are two main options for wetting and drying code (wd):
1345(a) an iterative limiter (il) and (b) a directional limiter (dl).
1346The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO
1347MS
1348which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative
1349limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
1350and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
1351by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.
1352
1353\begin{listing}
1354  \nlst{namwad}
1355  \caption{\forcode{&namwad}}
1356  \label{lst:namwad}
1357\end{listing}
1358
1359The following terminology is used. The depth of the topography (positive downwards)
1360at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code.
1361The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign
1362conventions used, the water depth, $h$, is the height of the free surface plus the depth of the
1363topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).
1364
1365Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
1366covered by water. They require the topography specified with a model
1367configuration to have negative depths at points where the land is higher than the
1368topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an
1369initial state with zero sea surface height elevation.
1370The user can choose to compute the vertical grid and heights in the model relative to
1371a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths
1372(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).
1373
1374Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.
1375$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies
1376with very steep slopes require larger values for normal choices of time-step. Surface fluxes
1377are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.
1378The fluxes are tappered down using a $\mathrm{tanh}$ weighting function
1379to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.
1380The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.
1381The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.
1382
1383Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration
1384and in ``realistic'' configurations covering parts of the north-west European shelf.
1385All these configurations have used pure sigma coordinates. It is expected that
1386the wetting and drying code will work in domains with more general s-coordinates provided
1387the coordinates are pure sigma in the region where wetting and drying actually occurs.
1388
1389The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.
1390The final sub-section covers some additional considerations that are relevant to both schemes.
1391
1392
1393%-----------------------------------------------------------------------------------------
1394%   Iterative limiters
1395%-----------------------------------------------------------------------------------------
1396\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})}
1397\label{subsec:DYN_wd_directional_limiter}
1398
1399The principal idea of the directional limiter is that
1400water 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}).
1401
1402All the changes associated with this option are made to the barotropic solver for the non-linear
1403free surface code within dynspg\_ts.
1404On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells
1405and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells
1406whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells
1407to be at least the minimum depth and hence is able to conserve mass / volume.
1408
1409The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).
1410If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the
1411flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets
1412\np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases
1413from 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.
1414
1415At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen
1416also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have
1417chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.
1418The rationale for setting the velocity to zero is that it is the momentum equations that are being solved
1419and the total momentum of the upstream cell (treating it as a finite volume) should be considered
1420to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its
1421treatment in the calculation of the flux of mass across the cell face.
1422
1423
1424\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic
1425timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
1426or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
1427fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
1428the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
1429to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the
1430baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.
1431
1432%-----------------------------------------------------------------------------------------
1433%   Iterative limiters
1434%-----------------------------------------------------------------------------------------
1435
1436\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})}
1437\label{subsec:DYN_wd_iterative_limiter}
1438
1439\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})}
1440\label{subsec:DYN_wd_il_spg_limiter}
1441
1442The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
1443or may become dry within the next time-step using an iterative method.
1444
1445The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:
1446
1447The continuity equation for the total water depth in a column
1448\begin{equation}
1449  \label{eq:DYN_wd_continuity}
1450  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
1451\end{equation}
1452can be written in discrete form  as
1453
1454\begin{align}
1455  \label{eq:DYN_wd_continuity_2}
1456  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
1457  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
1458  &= \mathrm{zzflx}_{i,j} .
1459\end{align}
1460
1461In the above $h$ is the depth of the water in the column at point $(i,j)$,
1462$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and
1463$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is
1464the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =
1465t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and
1466$\mathrm{zzflx}$ is the sum of the fluxes through all the faces.
1467
1468The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell
1469(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly
1470
1471\begin{equation}
1472  \label{eq:DYN_wd_zzflx_p_n_1}
1473  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
1474\end{equation}
1475
1476The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until
1477none of the cells will ``dry out''. To be precise the fluxes are limited until none of the
1478cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.
1479
1480Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and
1481$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients,
1482$\mathrm{zcoef}_{i,j}^{(m)}$ such that:
1483
1484\begin{equation}
1485  \label{eq:DYN_wd_continuity_coef}
1486  \begin{split}
1487    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
1488    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
1489  \end{split}
1490\end{equation}
1491
1492where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around
1493cells that would otherwise dry.
1494
1495The iteration is initialised by setting
1496
1497\begin{equation}
1498  \label{eq:DYN_wd_zzflx_initial}
1499  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
1500\end{equation}
1501
1502The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
1503cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
1504times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this
1505condition is
1506
1507\begin{equation}
1508  \label{eq:DYN_wd_continuity_if}
1509  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} ) .
1510\end{equation}
1511
1512Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum
1513outward flux that can be allowed and still maintain the minimum wet depth:
1514
1515\begin{equation}
1516  \label{eq:DYN_wd_max_flux}
1517  \begin{split}
1518    \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{]} \\
1519    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
1520  \end{split}
1521\end{equation}
1522
1523Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is
1524this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an
1525expression for the coefficient needed to multiply the outward flux at this cell in order
1526to avoid drying.
1527
1528\begin{equation}
1529  \label{eq:DYN_wd_continuity_nxtcoef}
1530  \begin{split}
1531    \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{]} \\
1532    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
1533  \end{split}
1534\end{equation}
1535
1536Only the outward flux components are altered but, of course, outward fluxes from one cell
1537are inward fluxes to adjacent cells and the balance in these cells may need subsequent
1538adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux
1539across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration
1540if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then
1541the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by
1542the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across
1543the faces of the cells do not ``compete'' (they do not over-write each other) and one
1544would expect the scheme to converge relatively quickly. The scheme is flux based so
1545conserves mass. It also conserves constant tracers for the same reason that the
1546directional limiter does.
1547
1548
1549%----------------------------------------------------------------------------------------
1550%      Surface pressure gradients
1551%----------------------------------------------------------------------------------------
1552\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})}
1553\label{subsec:DYN_wd_il_spg}
1554
1555At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
1556topography is sloping at these points the sea-surface will have a similar slope and there
1557will hence be very large horizontal pressure gradients at these points. The WAD modifies
1558the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such
1559points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie
1560between $0$ and $1$.
1561
1562We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated
1563at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the
1564neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
1565variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
1566column.  The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}.
1567
1568%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1569\begin{figure}[!ht]
1570  \centering
1571  \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg}
1572  \caption[Combinations controlling the limiting of the horizontal pressure gradient in
1573  wetting and drying regimes]{
1574    Three possible combinations of the logical variables controlling the
1575    limiting of the horizontal pressure gradient in wetting and drying regimes}
1576  \label{fig:DYN_WAD_dynhpg}
1577\end{figure}
1578%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1579
1580The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at
1581both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and
1582the minimum height of the sea surface at the two points is greater than the maximum height
1583of the topography at the two points:
1584
1585\begin{equation}
1586  \label{eq:DYN_ll_tmp1}
1587  \begin{split}
1588    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
1589                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
1590                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
1591                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
1592                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
1593  \end{split}
1594\end{equation}
1595
1596The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height
1597of the sea surface at the two points is greater than the maximum height of the topography
1598at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$
1599
1600\begin{equation}
1601  \label{eq:DYN_ll_tmp2}
1602  \begin{split}
1603    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
1604    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
1605    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
1606  \end{split}
1607\end{equation}
1608
1609If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$
1610point is unmodified. If both logicals are false zhpi is set to zero.
1611
1612If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure
1613gradient is multiplied through by zcpx which is the absolute value of the difference in
1614the water depths at the two points divided by the difference in the surface heights at the
1615two points. Thus the sign of the sea surface height gradient is retained but the magnitude
1616of the pressure force is determined by the difference in water depths rather than the
1617difference in surface height between the two points. Note that dividing by the difference
1618between the sea surface heights can be problematic if the heights approach parity. An
1619additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such
1620conditions.
1621
1622\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})}
1623\label{subsec:DYN_WAD_additional}
1624
1625In the very shallow water where wetting and drying occurs the parametrisation of
1626bottom drag is clearly very important. In order to promote stability
1627it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.
1628
1629Suitable specifcation of the surface heat flux in wetting and drying domains in forced and
1630coupled simulations needs further consideration. In order to prevent freezing or boiling
1631in uncoupled integrations the net surface heat fluxes need to be appropriately limited.
1632
1633%----------------------------------------------------------------------------------------
1634%      The WAD test cases
1635%----------------------------------------------------------------------------------------
1636\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})}
1637\label{subsec:DYN_WAD_test_cases}
1638
1639See the WAD tests MY\_DOC documention for details of the WAD test cases.
1640
1641
1642
1643% ================================================================
1644% Time evolution term
1645% ================================================================
1646\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}
1647\label{sec:DYN_nxt}
1648
1649%----------------------------------------------namdom----------------------------------------------------
1650
1651%-------------------------------------------------------------------------------------------------------------
1652
1653Options are defined through the \nam{dom}{dom} namelist variables.
1654The general framework for dynamics time stepping is a leap-frog scheme,
1655\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}).
1656The scheme is applied to the velocity, except when
1657using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})
1658in the variable volume case (\texttt{vvl?} defined),
1659where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})
1660
1661$\bullet$ vector invariant form or linear free surface
1662(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} not defined):
1663\[
1664  % \label{eq:DYN_nxt_vec}
1665  \left\{
1666    \begin{aligned}
1667      &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\
1668      &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
1669    \end{aligned}
1670  \right.
1671\]
1672
1673$\bullet$ flux form and nonlinear free surface
1674(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} ; \texttt{vvl?} defined):
1675\[
1676  % \label{eq:DYN_nxt_flux}
1677  \left\{
1678    \begin{aligned}
1679      &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\
1680      &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
1681      +\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]
1682    \end{aligned}
1683  \right.
1684\]
1685where RHS is the right hand side of the momentum equation,
1686the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.
1687$\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter).
1688Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}.
1689In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for
1690the momentum equations.
1691
1692Note that with the filtered free surface,
1693the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module,
1694and only array swapping and Asselin filtering is done in \mdl{dynnxt}.
1695
1696% ================================================================
1697\onlyinsubfile{\bibliography{../main/bibliography}}
1698
1699\onlyinsubfile{\printindex}
1700
1701\end{document}
Note: See TracBrowser for help on using the repository browser.