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

Last change on this file since 11690 was 11690, checked in by nicolasmartin, 12 months ago

Update chapters according to previous commit

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