source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_DYN.tex @ 10419

Last change on this file since 10419 was 10419, checked in by smasson, 21 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

File size: 73.7 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4% ================================================================
5% Chapter ——— Ocean Dynamics (DYN)
6% ================================================================
7\chapter{Ocean Dynamics (DYN)}
8\label{chap:DYN}
9
10\minitoc
11
12Using the representation described in \autoref{chap:DOM},
13several semi-discrete space forms of the dynamical equations are available depending on
14the vertical coordinate used and on the conservation properties of the vorticity term.
15In all the equations presented here, the masking has been omitted for simplicity.
16One must be aware that all the quantities are masked fields and
17that each time an average or difference operator is used, the resulting field is multiplied by a mask.
18
19The prognostic ocean dynamics equation can be summarized as follows:
20\[
21  \text{NXT} = \dbinom  {\text{VOR} + \text{KEG} + \text {ZAD} }
22  {\text{COR} + \text{ADV}                       }
23  + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
24\]
25NXT stands for next, referring to the time-stepping.
26The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that
27are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and
28a vertical advection part (ZAD) in the vector invariant formulation,
29or a Coriolis and advection part (COR+ADV) in the flux formulation.
30The terms following these are the pressure gradient contributions
31(HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient);
32and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF),
33which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules.
34The vertical diffusion term includes the surface and bottom stresses.
35The external forcings and parameterisations require complex inputs
36(surface wind stress calculation using bulk formulae, estimation of mixing coefficients)
37that are carried out in modules SBC, LDF and ZDF and are described in
38\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.
39
40In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence,
41curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module).
42
43The different options available to the user are managed by namelist variables.
44For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx},
45where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
46If a CPP key is used for this term its name is \key{ttt}.
47The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory,
48and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine.
49
50The user has the option of extracting and outputting each tendency term from the 3D momentum equations
51(\key{trddyn} defined), as described in \autoref{chap:MISC}.
52Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \key{trdvor} is defined)
53can be derived from the 3D terms.
54%%%
55\gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does
56MISC correspond to "extracting tendency terms" or "vorticity balance"?}
57
58% ================================================================
59% Sea Surface Height evolution & Diagnostics variables
60% ================================================================
61\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)}
62\label{sec:DYN_divcur_wzv}
63
64%--------------------------------------------------------------------------------------------------------------
65%           Horizontal divergence and relative vorticity
66%--------------------------------------------------------------------------------------------------------------
67\subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
68\label{subsec:DYN_divcur}
69
70The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows:
71\begin{equation}
72  \label{eq:divcur_cur}
73  \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right]
74      -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right)
75\end{equation} 
76
77The horizontal divergence is defined at a $T$-point.
78It is given by:
79\[
80  % \label{eq:divcur_div}
81  \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
82  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right]
83      +\delta_j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right)
84\]
85
86Note that although the vorticity has the same discrete expression in $z$- and $s$-coordinates,
87its physical meaning is not identical.
88$\zeta$ is a pseudo vorticity along $s$-surfaces
89(only pseudo because $(u,v)$ are still defined along geopotential surfaces,
90but are not necessarily defined at the same depth).
91
92The vorticity and divergence at the \textit{before} step are used in the computation of
93the horizontal diffusion of momentum.
94Note that because they have been calculated prior to the Asselin filtering of the \textit{before} velocities,
95the \textit{before} vorticity and divergence arrays must be included in the restart file to
96ensure perfect restartability.
97The vorticity and divergence at the \textit{now} time step are used for the computation of
98the nonlinear advection and of the vertical velocity respectively.
99
100%--------------------------------------------------------------------------------------------------------------
101%           Sea Surface Height evolution
102%--------------------------------------------------------------------------------------------------------------
103\subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
104\label{subsec:DYN_sshwzv}
105
106The sea surface height is given by:
107\begin{equation}
108  \label{eq:dynspg_ssh}
109  \begin{aligned}
110    \frac{\partial \eta }{\partial t}
111    &\equiv    \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{  \delta_i \left[ {e_{2u}\,e_{3u}\;u} \right]
112        +\delta_j \left[ {e_{1v}\,e_{3v}\;v} \right]  \right\} }
113    -    \frac{\textit{emp}}{\rho_w }   \\
114    &\equiv    \sum\limits_k {\chi \ e_{3t}}  -  \frac{\textit{emp}}{\rho_w }
115  \end{aligned}
116\end{equation}
117where \textit{emp} is the surface freshwater budget (evaporation minus precipitation),
118expressed in Kg/m$^2$/s (which is equal to mm/s),
119and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation).
120If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then
121\textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.
122The sea-surface height is evaluated using exactly the same time stepping scheme as
123the tracer equation \autoref{eq:tra_nxt}:
124a leapfrog scheme in combination with an Asselin time filter,
125$i.e.$ the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity).
126This is of paramount importance.
127Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
128the sea surface height equation otherwise tracer content will not be conserved
129\citep{Griffies_al_MWR01, Leclair_Madec_OM09}.
130
131The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
132taking into account the change of the thickness of the levels:
133\begin{equation}
134  \label{eq:wzv}
135  \left\{
136    \begin{aligned}
137      &\left. w \right|_{k_b-1/2} \quad= 0    \qquad \text{where } k_b \text{ is the level just above the sea floor }   \\
138      &\left. w \right|_{k+1/2}     = \left. w \right|_{k-1/2}  +  \left. e_{3t} \right|_{k}\;  \left. \chi \right|_k
139      - \frac{1} {2 \rdt} \left\left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right)
140    \end{aligned}
141  \right.
142\end{equation}
143
144In the case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$,
145as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses,
146re-orientated downward.
147\gmcomment{not sure of this...  to be modified with the change in emp setting}
148In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears.
149The upper boundary condition applies at a fixed level $z=0$.
150The top vertical velocity is thus equal to the divergence of the barotropic transport
151($i.e.$ the first term in the right-hand-side of \autoref{eq:dynspg_ssh}).
152
153Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
154its physical meaning is not the same:
155in the second case, $w$ is the velocity normal to the $s$-surfaces.
156Note also that the $k$-axis is re-orientated downwards in the \textsc{fortran} code compared to
157the indexing used in the semi-discrete equations such as \autoref{eq:wzv}
158(see \autoref{subsec:DOM_Num_Index_vertical}).
159
160
161% ================================================================
162% Coriolis and Advection terms: vector invariant form
163% ================================================================
164\section{Coriolis and advection: vector invariant form}
165\label{sec:DYN_adv_cor_vect}
166%-----------------------------------------nam_dynadv----------------------------------------------------
167
168\nlst{namdyn_adv} 
169%-------------------------------------------------------------------------------------------------------------
170
171The vector invariant form of the momentum equations is the one most often used in
172applications of the \NEMO ocean model.
173The flux form option (see next section) has been present since version $2$.
174Options are defined through the \ngn{namdyn\_adv} namelist variables Coriolis and
175momentum advection terms are evaluated using a leapfrog scheme,
176$i.e.$ the velocity appearing in these expressions is centred in time (\textit{now} velocity).
177At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following
178\autoref{chap:LBC}.
179
180% -------------------------------------------------------------------------------------------------------------
181%        Vorticity term
182% -------------------------------------------------------------------------------------------------------------
183\subsection{Vorticity term (\protect\mdl{dynvor})}
184\label{subsec:DYN_vor}
185%------------------------------------------nam_dynvor----------------------------------------------------
186
187\nlst{namdyn_vor} 
188%-------------------------------------------------------------------------------------------------------------
189
190Options are defined through the \ngn{namdyn\_vor} namelist variables.
191Four discretisations of the vorticity term (\np{ln\_dynvor\_xxx}\forcode{ = .true.}) are available:
192conserving potential enstrophy of horizontally non-divergent flow (ENS scheme);
193conserving horizontal kinetic energy (ENE scheme);
194conserving potential enstrophy for the relative vorticity term and
195horizontal kinetic energy for the planetary vorticity term (MIX scheme);
196or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy
197(EEN scheme) (see \autoref{subsec:C_vorEEN}).
198In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of
199vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}).
200The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module.
201
202%-------------------------------------------------------------
203%                 enstrophy conserving scheme
204%-------------------------------------------------------------
205\subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
206\label{subsec:DYN_vor_ens}
207
208In the enstrophy conserving case (ENS scheme),
209the discrete formulation of the vorticity term provides a global conservation of the enstrophy
210($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow ($i.e.$ $\chi$=$0$),
211but does not conserve the total kinetic energy.
212It is given by:
213\begin{equation}
214  \label{eq:dynvor_ens}
215  \left\{
216    \begin{aligned}
217      {+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i}
218      & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2}    \\
219      {- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}
220      & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j}
221    \end{aligned}
222  \right.
223\end{equation} 
224
225%-------------------------------------------------------------
226%                 energy conserving scheme
227%-------------------------------------------------------------
228\subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
229\label{subsec:DYN_vor_ene}
230
231The kinetic energy conserving scheme (ENE scheme) conserves the global kinetic energy but not the global enstrophy.
232It is given by:
233\begin{equation}
234  \label{eq:dynvor_ene}
235  \left\{
236    \begin{aligned}
237      {+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
238            \;  \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
239      {- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
240            \;  \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} }
241    \end{aligned}
242  \right.
243\end{equation} 
244
245%-------------------------------------------------------------
246%                 mix energy/enstrophy conserving scheme
247%-------------------------------------------------------------
248\subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) }
249\label{subsec:DYN_vor_mix}
250
251For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used.
252It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term,
253and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term.
254\[
255  % \label{eq:dynvor_mix}
256  \left\{ {
257      \begin{aligned}
258        {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i}
259          \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} }
260          \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
261              \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
262        {-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j
263          \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} }
264          \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
265              \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill
266      \end{aligned}
267    } \right.
268\]
269
270%-------------------------------------------------------------
271%                 energy and enstrophy conserving scheme
272%-------------------------------------------------------------
273\subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) }
274\label{subsec:DYN_vor_een}
275
276In both the ENS and ENE schemes,
277it is apparent that the combination of $i$ and $j$ averages of the velocity allows for
278the presence of grid point oscillation structures that will be invisible to the operator.
279These structures are \textit{computational modes} that will be at least partly damped by
280the momentum diffusion operator ($i.e.$ the subgrid-scale advection), but not by the resolved advection term.
281The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field.
282Such noise would result in more noise in the vertical velocity field, an undesirable feature.
283This is a well-known characteristic of $C$-grid discretization where
284$u$ and $v$ are located at different grid points,
285a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid.
286\gmcomment{ To circumvent this, Adcroft (ADD REF HERE)
287Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
288
289A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}.
290The idea is to get rid of the double averaging by considering triad combinations of vorticity.
291It is noteworthy that this solution is conceptually quite similar to the one proposed by
292\citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
293
294The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified
295for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.
296First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
297\[
298  % \label{eq:pot_vor}
299  q  = \frac{\zeta +f} {e_{3f} }
300\]
301where the relative vorticity is defined by (\autoref{eq:divcur_cur}),
302the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
303\begin{equation}
304  \label{eq:een_e3f}
305  e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
306\end{equation}
307
308%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
309\begin{figure}[!ht]
310  \begin{center}
311    \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad}
312    \caption{
313      \protect\label{fig:DYN_een_triad}
314      Triads used in the energy and enstrophy conserving scheme (een) for
315      $u$-component (upper panel) and $v$-component (lower panel).
316    }
317  \end{center}
318\end{figure}
319% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
320
321A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.
322It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks
323(\np{nn\_een\_e3f}\forcode{ = 1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{ = .true.}).
324The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
325extends by continuity the value of $e_{3f}$ into the land areas.
326This case introduces a sub-grid-scale topography at f-points
327(with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry)
328that tends to reinforce the topostrophy of the flow
329($i.e.$ the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.
330
331Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
332the following triad combinations of the neighbouring potential vorticities defined at f-points
333(\autoref{fig:DYN_een_triad}):
334\begin{equation}
335  \label{eq:Q_triads}
336  _i^j \mathbb{Q}^{i_p}_{j_p}
337  = \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)
338\end{equation}
339where 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$.
340
341Finally, the vorticity terms are represented as:
342\begin{equation}
343  \label{eq:dynvor_een}
344  \left\{ {
345      \begin{aligned}
346        +q\,e_3 \, v    &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}}
347        {^{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}   \\
348        - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}}
349        {^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}   \\
350      \end{aligned}
351    } \right.
352\end{equation} 
353
354This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
355It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow
356($i.e.$ $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).
357Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
358the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.
359Furthermore, used in combination with a partial steps representation of bottom topography,
360it improves the interaction between current and topography,
361leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}.
362
363%--------------------------------------------------------------------------------------------------------------
364%           Kinetic Energy Gradient term
365%--------------------------------------------------------------------------------------------------------------
366\subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})}
367\label{subsec:DYN_keg}
368
369As demonstrated in \autoref{apdx:C},
370there is a single discrete formulation of the kinetic energy gradient term that,
371together with the formulation chosen for the vertical advection (see below),
372conserves the total kinetic energy:
373\[
374  % \label{eq:dynkeg}
375  \left\{
376    \begin{aligned}
377      -\frac{1}{2 \; e_{1u} }  & \ \delta_{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\
378      -\frac{1}{2 \; e_{2v} }  & \ \delta_{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]
379    \end{aligned}
380  \right.
381\]
382
383%--------------------------------------------------------------------------------------------------------------
384%           Vertical advection term
385%--------------------------------------------------------------------------------------------------------------
386\subsection{Vertical advection term (\protect\mdl{dynzad}) }
387\label{subsec:DYN_zad}
388
389The discrete formulation of the vertical advection, t
390ogether with the formulation chosen for the gradient of kinetic energy (KE) term,
391conserves the total kinetic energy.
392Indeed, the change of KE due to the vertical advection is exactly balanced by
393the change of KE due to the gradient of KE (see \autoref{apdx:C}).
394\[
395  % \label{eq:dynzad}
396  \left\{
397    \begin{aligned}
398      -\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}  \\
399      -\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}
400    \end{aligned}
401  \right.
402\]
403When \np{ln\_dynzad\_zts}\forcode{ = .true.},
404a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term.
405This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
406Note that in this case,
407a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability,
408an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}).
409
410
411% ================================================================
412% Coriolis and Advection : flux form
413% ================================================================
414\section{Coriolis and advection: flux form}
415\label{sec:DYN_adv_cor_flux}
416%------------------------------------------nam_dynadv----------------------------------------------------
417
418\nlst{namdyn_adv} 
419%-------------------------------------------------------------------------------------------------------------
420
421Options are defined through the \ngn{namdyn\_adv} namelist variables.
422In the flux form (as in the vector invariant form),
423the Coriolis and momentum advection terms are evaluated using a leapfrog scheme,
424$i.e.$ the velocity appearing in their expressions is centred in time (\textit{now} velocity).
425At the lateral boundaries either free slip,
426no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}.
427
428
429%--------------------------------------------------------------------------------------------------------------
430%           Coriolis plus curvature metric terms
431%--------------------------------------------------------------------------------------------------------------
432\subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) }
433\label{subsec:DYN_cor_flux}
434
435In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term.
436This altered Coriolis parameter is thus discretised at $f$-points.
437It is given by:
438\begin{multline*}
439  % \label{eq:dyncor_metric}
440  f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right\\
441  \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right]
442      -  \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right]  }  \ \right)
443\end{multline*} 
444
445Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to
446compute the product of the Coriolis parameter and the vorticity.
447However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date.
448This term is evaluated using a leapfrog scheme, $i.e.$ the velocity is centred in time (\textit{now} velocity).
449
450%--------------------------------------------------------------------------------------------------------------
451%           Flux form Advection term
452%--------------------------------------------------------------------------------------------------------------
453\subsection{Flux form advection term (\protect\mdl{dynadv}) }
454\label{subsec:DYN_adv_flux}
455
456The discrete expression of the advection term is given by:
457\[
458  % \label{eq:dynadv}
459  \left\{
460    \begin{aligned}
461      \frac{1}{e_{1u}\,e_{2u}\,e_{3u}}
462      \left(      \delta_{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i       }  \ u_t      \right]
463        + \delta_{j       } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2\ u_f      \right] \right\ \;   \\
464      \left.   + \delta_{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2\ u_{uw} \right] \right)   \\
465      \\
466      \frac{1}{e_{1v}\,e_{2v}\,e_{3v}}
467      \left(     \delta_{i       } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f       \right]
468        + \delta_{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i       } \ v_t       \right] \right\ \, \, \\
469      \left+ \delta_{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw}  \right] \right) \\
470    \end{aligned}
471  \right.
472\]
473
474Two advection schemes are available:
475a $2^{nd}$ order centered finite difference scheme, CEN2,
476or a $3^{rd}$ order upstream biased scheme, UBS.
477The latter is described in \citet{Shchepetkin_McWilliams_OM05}.
478The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.
479In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
480$u$ and $v$ at the centre of each face of $u$- and $v$-cells, $i.e.$ at the $T$-, $f$-,
481and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$.
482
483%-------------------------------------------------------------
484%                 2nd order centred scheme
485%-------------------------------------------------------------
486\subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}
487\label{subsec:DYN_adv_cen2}
488
489In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points:
490\begin{equation}
491  \label{eq:dynadv_cen2}
492  \left\{
493    \begin{aligned}
494      u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\
495      v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j    \quad &  v_{vw}^{cen2} &=\overline v ^{k+1/2}  \\
496    \end{aligned}
497  \right.
498\end{equation} 
499
500The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive ($i.e.$ it may create false extrema).
501It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to
502produce a sensible solution.
503The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter,
504so $u$ and $v$ are the \emph{now} velocities.
505
506%-------------------------------------------------------------
507%                 UBS scheme
508%-------------------------------------------------------------
509\subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}
510\label{subsec:DYN_adv_ubs}
511
512The UBS advection scheme is an upstream biased third order scheme based on
513an upstream-biased parabolic interpolation.
514For example, the evaluation of $u_T^{ubs} $ is done as follows:
515\begin{equation}
516  \label{eq:dynadv_ubs}
517  u_T^{ubs} =\overline u ^i-\;\frac{1}{6}
518  \begin{cases}
519    u"_{i-1/2}&   \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  \geqslant 0$ }    \\
520    u"_{i+1/2}&   \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  < 0$ }
521  \end{cases}
522\end{equation}
523where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$.
524This results in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error
525\citep{Shchepetkin_McWilliams_OM05}.
526The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}.
527It is a relatively good compromise between accuracy and smoothness.
528It is not a \emph{positive} scheme, meaning that false extrema are permitted.
529But the amplitudes of the false extrema are significantly reduced over those in the centred second order method.
530As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum
531($i.e.$ \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ = .false.}),
532and it is recommended to do so.
533
534The UBS scheme is not used in all directions.
535In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and
536$u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used.
537UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the
538sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }
539
540For stability reasons, the first term in (\autoref{eq:dynadv_ubs}),
541which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time),
542while the second term, which is the diffusion part of the scheme,
543is evaluated using the \textit{before} velocity (forward in time).
544This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
545
546Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
547one coefficient.
548Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
549This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
550Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
551
552Note also that in the current version of \mdl{dynadv\_ubs},
553there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS.
554This is an error and should be suppressed soon.
555%%%
556\gmcomment{action :  this have to be done}
557%%%
558
559% ================================================================
560%           Hydrostatic pressure gradient term
561% ================================================================
562\section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
563\label{sec:DYN_hpg}
564%------------------------------------------nam_dynhpg---------------------------------------------------
565
566\nlst{namdyn_hpg} 
567%-------------------------------------------------------------------------------------------------------------
568
569Options are defined through the \ngn{namdyn\_hpg} namelist variables.
570The key distinction between the different algorithms used for
571the hydrostatic pressure gradient is the vertical coordinate used,
572since HPG is a \emph{horizontal} pressure gradient, $i.e.$ computed along geopotential surfaces.
573As a result, any tilt of the surface of the computational levels will require a specific treatment to
574compute the hydrostatic pressure gradient.
575
576The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme,
577$i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$),
578or a semi-implcit scheme.
579At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied.
580
581%--------------------------------------------------------------------------------------------------------------
582%           z-coordinate with full step
583%--------------------------------------------------------------------------------------------------------------
584\subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}
585\label{subsec:DYN_hpg_zco}
586
587The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface.
588However, the pressure is large at great depth while its horizontal gradient is several orders of magnitude smaller.
589This may lead to large truncation errors in the pressure gradient terms.
590Thus, the two horizontal components of the hydrostatic pressure gradient are computed directly as follows:
591
592for $k=km$ (surface layer, $jk=1$ in the code)
593\begin{equation}
594  \label{eq:dynhpg_zco_surf}
595  \left\{
596    \begin{aligned}
597      \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k=km}
598      &= \frac{1}{2} g \   \left. \delta_{i+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
599      \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k=km}
600      &= \frac{1}{2} g \   \left. \delta_{j+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
601    \end{aligned}
602  \right.
603\end{equation} 
604
605for $1<k<km$ (interior layer)
606\begin{equation}
607  \label{eq:dynhpg_zco}
608  \left\{
609    \begin{aligned}
610      \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k}
611      &=             \left. \delta_{i+1/2} \left[  p^h          \right] \right|_{k-1}
612      +    \frac{1}{2}\;g\;   \left. \delta_{i+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
613      \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k}
614      &=                \left. \delta_{j+1/2} \left[  p^h          \right] \right|_{k-1}
615      +    \frac{1}{2}\;g\;   \left. \delta_{j+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
616    \end{aligned}
617  \right.
618\end{equation} 
619
620Note that the $1/2$ factor in (\autoref{eq:dynhpg_zco_surf}) is adequate because of the definition of $e_{3w}$ as
621the vertical derivative of the scale factor at the surface level ($z=0$).
622Note also that in case of variable volume level (\key{vvl} defined),
623the surface pressure gradient is included in \autoref{eq:dynhpg_zco_surf} and
624\autoref{eq:dynhpg_zco} through the space and time variations of the vertical scale factor $e_{3w}$.
625
626%--------------------------------------------------------------------------------------------------------------
627%           z-coordinate with partial step
628%--------------------------------------------------------------------------------------------------------------
629\subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}
630\label{subsec:DYN_hpg_zps}
631
632With partial bottom cells, tracers in horizontally adjacent cells generally live at different depths.
633Before taking horizontal gradients between these tracer points,
634a linear interpolation is used to approximate the deeper tracer as if
635it actually lived at the depth of the shallower tracer point.
636
637Apart from this modification,
638the horizontal hydrostatic pressure gradient evaluated in the $z$-coordinate with partial step is exactly as in
639the pure $z$-coordinate case.
640As explained in detail in section \autoref{sec:TRA_zpshde},
641the nonlinearity of pressure effects in the equation of state is such that
642it is better to interpolate temperature and salinity vertically before computing the density.
643Horizontal gradients of temperature and salinity are needed for the TRA modules,
644which is the reason why the horizontal gradients of density at the deepest model level are computed in
645module \mdl{zpsdhe} located in the TRA directory and described in \autoref{sec:TRA_zpshde}.
646
647%--------------------------------------------------------------------------------------------------------------
648%           s- and s-z-coordinates
649%--------------------------------------------------------------------------------------------------------------
650\subsection{$S$- and $Z$-$S$-coordinates}
651\label{subsec:DYN_hpg_sco}
652
653Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers
654($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}).
655A number of different pressure gradient options are coded but the ROMS-like,
656density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
657
658$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
659\begin{equation}
660  \label{eq:dynhpg_sco}
661  \left\{
662    \begin{aligned}
663      - \frac{1}                 {\rho_o \, e_{1u}} \;   \delta_{i+1/2} \left[  p^h  \right]
664      + \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta_{i+1/2} \left[  z_t   \right]    \\
665      - \frac{1}                 {\rho_o \, e_{2v}} \;   \delta_{j+1/2} \left[  p^h  \right]
666      + \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta_{j+1/2} \left[  z_t   \right]    \\
667    \end{aligned}
668  \right.
669\end{equation} 
670
671Where the first term is the pressure gradient along coordinates,
672computed as in \autoref{eq:dynhpg_zco_surf} - \autoref{eq:dynhpg_zco},
673and $z_T$ is the depth of the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
674($e_{3w}$).
675 
676$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}\forcode{ = .true.}).
677This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}\forcode{ = .true.}).
678
679$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.})
680
681$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} 
682(\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development)
683
684Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated
685(\key{vvl}) because in that case, even with a flat bottom,
686the coordinate surfaces are not horizontal but follow the free surface \citep{Levier2007}.
687The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as
688an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active.
689The pressure Jacobian scheme uses a constrained cubic spline to
690reconstruct the density profile across the water column.
691This method maintains the monotonicity between the density nodes.
692The pressure can be calculated by analytical integration of the density profile and
693a pressure Jacobian method is used to solve the horizontal pressure gradient.
694This method can provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.
695
696\subsection{Ice shelf cavity}
697\label{subsec:DYN_hpg_isf}
698Beneath an ice shelf, the total pressure gradient is the sum of the pressure gradient due to the ice shelf load and
699the pressure gradient due to the ocean load.
700If cavity opened (\np{ln\_isfcav}\forcode{ = .true.}) these 2 terms can be calculated by
701setting \np{ln\_dynhpg\_isf}\forcode{ = .true.}.
702No other scheme are working with the ice shelf.\\
703
704$\bullet$ The main hypothesis to compute the ice shelf load is that the ice shelf is in an isostatic equilibrium.
705The top pressure is computed integrating from surface to the base of the ice shelf a reference density profile
706(prescribed as density of a water at 34.4 PSU and -1.9\degC) and
707corresponds to the water replaced by the ice shelf.
708This top pressure is constant over time.
709A detailed description of this method is described in \citet{Losch2008}.\\
710
711$\bullet$ The ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in
712\autoref{subsec:DYN_hpg_sco}.
713
714%--------------------------------------------------------------------------------------------------------------
715%           Time-scheme
716%--------------------------------------------------------------------------------------------------------------
717\subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})}
718\label{subsec:DYN_hpg_imp}
719
720The default time differencing scheme used for the horizontal pressure gradient is a leapfrog scheme and
721therefore the density used in all discrete expressions given above is the  \textit{now} density,
722computed from the \textit{now} temperature and salinity.
723In some specific cases
724(usually high resolution simulations over an ocean domain which includes weakly stratified regions)
725the physical phenomenon that controls the time-step is internal gravity waves (IGWs).
726A semi-implicit scheme for doubling the stability limit associated with IGWs can be used
727\citep{Brown_Campana_MWR78, Maltrud1998}.
728It involves the evaluation of the hydrostatic pressure gradient as
729an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$
730($i.e.$  \textit{before}\textit{now} and  \textit{after} time-steps),
731rather than at the central time level $t$ only, as in the standard leapfrog scheme.
732
733$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}):
734
735\begin{equation}
736  \label{eq:dynhpg_lf}
737  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
738  -\frac{1}{\rho_o \,e_{1u} }\delta_{i+1/2} \left[ {p_h^t } \right]
739\end{equation}
740
741$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}\forcode{ = .true.}):
742\begin{equation}
743  \label{eq:dynhpg_imp}
744  \frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
745  -\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]
746\end{equation}
747
748The semi-implicit time scheme \autoref{eq:dynhpg_imp} is made possible without
749significant additional computation since the density can be updated to time level $t+\rdt$ before
750computing the horizontal hydrostatic pressure gradient.
751It can be easily shown that the stability limit associated with the hydrostatic pressure gradient doubles using
752\autoref{eq:dynhpg_imp} compared to that using the standard leapfrog scheme \autoref{eq:dynhpg_lf}.
753Note that \autoref{eq:dynhpg_imp} is equivalent to applying a time filter to the pressure gradient to
754eliminate high frequency IGWs.
755Obviously, when using \autoref{eq:dynhpg_imp},
756the doubling of the time-step is achievable only if no other factors control the time-step,
757such as the stability limits associated with advection or diffusion.
758
759In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}\forcode{ = .true.}.
760In this case, we choose to apply the time filter to temperature and salinity used in the equation of state,
761instead of applying it to the hydrostatic pressure or to the density,
762so that no additional storage array has to be defined.
763The density used to compute the hydrostatic pressure gradient (whatever the formulation) is evaluated as follows:
764\[
765  % \label{eq:rho_flt}
766  \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
767  \quad    \text{with}  \quad
768  \widetilde{X} = 1 / 4 \left(  X^{t+\rdt} +2 \,X^t + X^{t-\rdt\right)
769\]
770
771Note that in the semi-implicit case, it is necessary to save the filtered density,
772an extra three-dimensional field, in the restart file to restart the model with exact reproducibility.
773This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter.
774
775% ================================================================
776% Surface Pressure Gradient
777% ================================================================
778\section{Surface pressure gradient (\protect\mdl{dynspg})}
779\label{sec:DYN_spg}
780%-----------------------------------------nam_dynspg----------------------------------------------------
781
782\nlst{namdyn_spg} 
783%------------------------------------------------------------------------------------------------------------
784
785Options are defined through the \ngn{namdyn\_spg} namelist variables.
786The surface pressure gradient term is related to the representation of the free surface (\autoref{sec:PE_hor_pg}).
787The main distinction is between the fixed volume case (linear free surface) and
788the variable volume case (nonlinear free surface, \key{vvl} is defined).
789In the linear free surface case (\autoref{subsec:PE_free_surface})
790the vertical scale factors $e_{3}$ are fixed in time,
791while they are time-dependent in the nonlinear case (\autoref{subsec:PE_free_surface}).
792With both linear and nonlinear free surface, external gravity waves are allowed in the equations,
793which imposes a very small time step when an explicit time stepping is used.
794Two methods are proposed to allow a longer time step for the three-dimensional equations:
795the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}),
796and the split-explicit free surface described below.
797The extra term introduced in the filtered method is calculated implicitly,
798so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
799
800
801The form of the surface pressure gradient term depends on how the user wants to
802handle the fast external gravity waves that are a solution of the analytical equation (\autoref{sec:PE_hor_pg}).
803Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx):
804an explicit formulation which requires a small time step;
805a filtered free surface formulation which allows a larger time step by
806adding a filtering term into the momentum equation;
807and a split-explicit free surface formulation, described below, which also allows a larger time step.
808
809The extra term introduced in the filtered method is calculated implicitly, so that a solver is used to compute it.
810As a consequence the update of the $next$ velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
811
812
813%--------------------------------------------------------------------------------------------------------------
814% Explicit free surface formulation
815%--------------------------------------------------------------------------------------------------------------
816\subsection{Explicit free surface (\protect\key{dynspg\_exp})}
817\label{subsec:DYN_spg_exp}
818
819In the explicit free surface formulation (\key{dynspg\_exp} defined),
820the model time step is chosen to be small enough to resolve the external gravity waves
821(typically a few tens of seconds).
822The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time),
823is thus simply given by :
824\begin{equation}
825  \label{eq:dynspg_exp}
826  \left\{
827    \begin{aligned}
828      - \frac{1}{e_{1u}\,\rho_o} \; \delta_{i+1/2} \left[  \,\rho \,\eta\,  \right]    \\
829      - \frac{1}{e_{2v}\,\rho_o} \; \delta_{j+1/2} \left[  \,\rho \,\eta\,  \right]
830    \end{aligned}
831  \right.
832\end{equation} 
833
834Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined),
835the surface pressure gradient is already included in the momentum tendency through
836the level thickness variation allowed in the computation of the hydrostatic pressure gradient.
837Thus, nothing is done in the \mdl{dynspg\_exp} module.
838
839%--------------------------------------------------------------------------------------------------------------
840% Split-explict free surface formulation
841%--------------------------------------------------------------------------------------------------------------
842\subsection{Split-explicit free surface (\protect\key{dynspg\_ts})}
843\label{subsec:DYN_spg_ts}
844%------------------------------------------namsplit-----------------------------------------------------------
845%
846%\nlst{namsplit}
847%-------------------------------------------------------------------------------------------------------------
848
849The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),
850also called the time-splitting formulation, follows the one proposed by \citet{Shchepetkin_McWilliams_OM05}.
851The general idea is to solve the free surface equation and the associated barotropic velocity equations with
852a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
853(\autoref{fig:DYN_dynspg_ts}).
854The size of the small time step, $\rdt_e$ (the external mode or barotropic time step) is provided through
855the \np{nn\_baro} namelist parameter as: $\rdt_e = \rdt / nn\_baro$.
856This parameter can be optionally defined automatically (\np{ln\_bt\_nn\_auto}\forcode{ = .true.}) considering that
857the stability of the barotropic system is essentially controled by external waves propagation.
858Maximum Courant number is in that case time independent, and easily computed online from the input bathymetry.
859Therefore, $\rdt_e$ is adjusted so that the Maximum allowed Courant number is smaller than \np{rn\_bt\_cmax}.
860
861%%%
862The barotropic mode solves the following equations:
863% \begin{subequations}
864%  \label{eq:BT}
865\begin{equation}
866  \label{eq:BT_dyn}
867  \frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}=
868  -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h}
869  -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}}
870\end{equation}
871\[
872  % \label{eq:BT_ssh}
873  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E
874\]
875% \end{subequations}
876where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes,
877surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency.
878The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress
879(see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration.
880Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
881detailed in \citet{Shchepetkin_McWilliams_OM05}.
882AB3-AM4 coefficients used in \NEMO follow the second-order accurate,
883"multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08}
884(see their figure 12, lower left).
885
886%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
887\begin{figure}[!t]
888  \begin{center}
889    \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts}
890    \caption{
891      \protect\label{fig:DYN_dynspg_ts}
892      Schematic of the split-explicit time stepping scheme for the external and internal modes.
893      Time increases to the right. In this particular exemple,
894      a boxcar averaging window over $nn\_baro$ barotropic time steps is used ($nn\_bt\_flt=1$) and $nn\_baro=5$.
895      Internal mode time steps (which are also the model time steps) are denoted by $t-\rdt$, $t$ and $t+\rdt$.
896      Variables with $k$ superscript refer to instantaneous barotropic variables,
897      $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary (red vertical bars) and
898      secondary weights (blue vertical bars).
899      The former are used to obtain time filtered quantities at $t+\rdt$ while
900      the latter are used to obtain time averaged transports to advect tracers.
901      a) Forward time integration: \protect\np{ln\_bt\_fw}\forcode{ = .true.},
902      \protect\np{ln\_bt\_av}\forcode{ = .true.}.
903      b) Centred time integration: \protect\np{ln\_bt\_fw}\forcode{ = .false.},
904      \protect\np{ln\_bt\_av}\forcode{ = .true.}.
905      c) Forward time integration with no time filtering (POM-like scheme):
906      \protect\np{ln\_bt\_fw}\forcode{ = .true.}, \protect\np{ln\_bt\_av}\forcode{ = .false.}.
907    }
908  \end{center}
909\end{figure}
910%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
911
912In the default case (\np{ln\_bt\_fw}\forcode{ = .true.}),
913the external mode is integrated between \textit{now} and \textit{after} baroclinic time-steps
914(\autoref{fig:DYN_dynspg_ts}a).
915To avoid aliasing of fast barotropic motions into three dimensional equations,
916time filtering is eventually applied on barotropic quantities (\np{ln\_bt\_av}\forcode{ = .true.}).
917In that case, the integration is extended slightly beyond \textit{after} time step to
918provide time filtered quantities.
919These are used for the subsequent initialization of the barotropic mode in the following baroclinic step.
920Since external mode equations written at baroclinic time steps finally follow a forward time stepping scheme,
921asselin filtering is not applied to barotropic quantities.\\
922Alternatively, one can choose to integrate barotropic equations starting from \textit{before} time step
923(\np{ln\_bt\_fw}\forcode{ = .false.}).
924Although more computationaly expensive ( \np{nn\_baro} additional iterations are indeed necessary),
925the baroclinic to barotropic forcing term given at \textit{now} time step become centred in
926the middle of the integration window.
927It can easily be shown that this property removes part of splitting errors between modes,
928which increases the overall numerical robustness.
929%references to Patrick Marsaleix' work here. Also work done by SHOM group.
930
931%%%
932
933As far as tracer conservation is concerned,
934barotropic velocities used to advect tracers must also be updated at \textit{now} time step.
935This implies to change the traditional order of computations in \NEMO:
936most of momentum trends (including the barotropic mode calculation) updated first, tracers' after.
937This \textit{de facto} makes semi-implicit hydrostatic pressure gradient
938(see section \autoref{subsec:DYN_hpg_imp})
939and time splitting not compatible.
940Advective barotropic velocities are obtained by using a secondary set of filtering weights,
941uniquely defined from the filter coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}).
942Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to
943obtain exact conservation.
944
945%%%
946
947One can eventually choose to feedback instantaneous values by not using any time filter
948(\np{ln\_bt\_av}\forcode{ = .false.}).
949In that case, external mode equations are continuous in time,
950$i.e.$ they are not re-initialized when starting a new sub-stepping sequence.
951This is the method used so far in the POM model, the stability being maintained by
952refreshing at (almost) each barotropic time step advection and horizontal diffusion terms.
953Since the latter terms have not been added in \NEMO for computational efficiency,
954removing time filtering is not recommended except for debugging purposes.
955This may be used for instance to appreciate the damping effect of the standard formulation on
956external gravity waves in idealized or weakly non-linear cases.
957Although the damping is lower than for the filtered free surface,
958it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.
959
960%>>>>>===============
961\gmcomment{               %%% copy from griffies Book
962
963\textbf{title: Time stepping the barotropic system }
964
965Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$.
966Hence, we can update the surface height and vertically integrated velocity with a leap-frog scheme using
967the small barotropic time step $\rdt$.
968We have
969
970\[
971  % \label{eq:DYN_spg_ts_eta}
972  \eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
973  = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right]
974\]
975\begin{multline*}
976  % \label{eq:DYN_spg_ts_u}
977  \textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}\\
978  = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})
979    - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right]
980\end{multline*}
981\
982
983In these equations, araised (b) denotes values of surface height and vertically integrated velocity updated with
984the barotropic time steps.
985The $\tau$ time label on $\eta^{(b)}$ and $U^{(b)}$ denotes the baroclinic time at which
986the vertically integrated forcing $\textbf{M}(\tau)$
987(note that this forcing includes the surface freshwater forcing),
988the tracer fields, the freshwater flux $\text{EMP}_w(\tau)$,
989and total depth of the ocean $H(\tau)$ are held for the duration of the barotropic time stepping over
990a single cycle.
991This is also the time that sets the barotropic time steps via
992\[
993  % \label{eq:DYN_spg_ts_t}
994  t_n=\tau+n\rdt
995\]
996with $n$ an integer.
997The density scaled surface pressure is evaluated via
998\[
999  % \label{eq:DYN_spg_ts_ps}
1000  p_s^{(b)}(\tau,t_{n}) =
1001  \begin{cases}
1002    g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o  &      \text{non-linear case} \\
1003    g \;\eta_s^{(b)}(\tau,t_{n})  &      \text{linear case}
1004  \end{cases}
1005\]
1006To get started, we assume the following initial conditions
1007\[
1008  % \label{eq:DYN_spg_ts_eta}
1009  \begin{split}
1010    \eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)}    \\
1011    \eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
1012  \end{split}
1013\]
1014with
1015\[
1016  % \label{eq:DYN_spg_ts_etaF}
1017  \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
1018\]
1019the time averaged surface height taken from the previous barotropic cycle.
1020Likewise,
1021\[
1022  % \label{eq:DYN_spg_ts_u}
1023  \textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)} \\ \\
1024  \textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}
1025\]
1026with
1027\[
1028  % \label{eq:DYN_spg_ts_u}
1029  \overline{\textbf{U}^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n})
1030\]
1031the time averaged vertically integrated transport.
1032Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration.
1033
1034Upon reaching $t_{n=N} = \tau + 2\rdt \tau$ ,
1035the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at
1036baroclinic time $\tau + \rdt \tau$ 
1037\[
1038  % \label{eq:DYN_spg_ts_u}
1039  \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})
1040\]
1041The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using
1042the following form
1043
1044\begin{equation}
1045  \label{eq:DYN_spg_ts_ssh}
1046  \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] 
1047\end{equation}
1048
1049The use of this "big-leap-frog" scheme for the surface height ensures compatibility between
1050the mass/volume budgets and the tracer budgets.
1051More discussion of this point is provided in Chapter 10 (see in particular Section 10.2).
1052 
1053In general, some form of time filter is needed to maintain integrity of the surface height field due to
1054the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
1055We have tried various forms of such filtering,
1056with the following method discussed in \cite{Griffies_al_MWR01} chosen due to
1057its stability and reasonably good maintenance of tracer conservation properties (see ??).
1058
1059\begin{equation}
1060  \label{eq:DYN_spg_ts_sshf}
1061  \eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)}
1062\end{equation}
1063Another approach tried was
1064
1065\[
1066  % \label{eq:DYN_spg_ts_sshf2}
1067  \eta^{F}(\tau-\Delta) = \eta(\tau)
1068  + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
1069    + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
1070\]
1071
1072which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$.
1073This isolation allows for an easy check that tracer conservation is exact when
1074eliminating tracer and surface height time filtering (see ?? for more complete discussion).
1075However, in the general case with a non-zero $\alpha$,
1076the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.
1077
1078}            %%end gm comment (copy of griffies book)
1079
1080%>>>>>===============
1081
1082
1083%--------------------------------------------------------------------------------------------------------------
1084% Filtered free surface formulation
1085%--------------------------------------------------------------------------------------------------------------
1086\subsection{Filtered free surface (\protect\key{dynspg\_flt})}
1087\label{subsec:DYN_spg_fltp}
1088
1089The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.
1090The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.
1091The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
1092
1093%% gm %%======>>>>   given here the discrete eqs provided to the solver
1094\gmcomment{               %%% copy from chap-model basics
1095  \[
1096    % \label{eq:spg_flt}
1097    \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
1098    - g \nabla \left( \tilde{\rho} \ \eta \right)
1099    - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)
1100  \]
1101  where $T_c$, is a parameter with dimensions of time which characterizes the force,
1102  $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density,
1103  and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
1104  non-linear and viscous terms in \autoref{eq:PE_dyn}.
1105}   %end gmcomment
1106
1107Note that in the linear free surface formulation (\key{vvl} not defined),
1108the ocean depth is time-independent and so is the matrix to be inverted.
1109It is computed once and for all and applies to all ocean time steps.
1110
1111% ================================================================
1112% Lateral diffusion term
1113% ================================================================
1114\section{Lateral diffusion term and operators (\protect\mdl{dynldf})}
1115\label{sec:DYN_ldf}
1116%------------------------------------------nam_dynldf----------------------------------------------------
1117
1118\nlst{namdyn_ldf} 
1119%-------------------------------------------------------------------------------------------------------------
1120
1121Options are defined through the \ngn{namdyn\_ldf} namelist variables.
1122The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators.
1123The coefficients may be constant or spatially variable;
1124the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}).
1125The lateral diffusion of momentum is evaluated using a forward scheme,
1126$i.e.$ the velocity appearing in its expression is the \textit{before} velocity in time,
1127except for the pure vertical component that appears when a tensor of rotation is used.
1128This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}).
1129
1130At the lateral boundaries either free slip,
1131no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}).
1132
1133\gmcomment{
1134  Hyperviscous operators are frequently used in the simulation of turbulent flows to
1135  control the dissipation of unresolved small scale features.
1136  Their primary role is to provide strong dissipation at the smallest scale supported by
1137  the grid while minimizing the impact on the larger scale features.
1138  Hyperviscous operators are thus designed to be more scale selective than the traditional,
1139  physically motivated Laplace operator.
1140  In finite difference methods,
1141  the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since
1142  its damping time ($i.e.$ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$
1143  (so that short waves damped more rapidelly than long ones),
1144  whereas the Laplace operator damping time scales only like $\lambda^{-2}$.
1145}
1146
1147% ================================================================
1148\subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})]
1149            {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})}
1150\label{subsec:DYN_ldf_lap}
1151
1152For lateral iso-level diffusion, the discrete operator is:
1153\begin{equation}
1154  \label{eq:dynldf_lap}
1155  \left\{
1156    \begin{aligned}
1157      D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
1158          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
1159        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
1160      D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
1161          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
1162        {A_f^{lm} \;e_{3f} \zeta } \right]
1163    \end{aligned}
1164  \right.
1165\end{equation} 
1166
1167As explained in \autoref{subsec:PE_ldf},
1168this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and
1169ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.
1170
1171%--------------------------------------------------------------------------------------------------------------
1172%           Rotated laplacian operator
1173%--------------------------------------------------------------------------------------------------------------
1174\subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})]
1175            {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})}
1176\label{subsec:DYN_ldf_iso}
1177
1178A rotation of the lateral momentum diffusion operator is needed in several cases:
1179for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}\forcode{ = .true.}) and
1180for either iso-neutral (\np{ln\_dynldf\_iso}\forcode{ = .true.}) or
1181geopotential (\np{ln\_dynldf\_hor}\forcode{ = .true.}) diffusion in the $s$-coordinate.
1182In the partial step case, coordinates are horizontal except at the deepest level and
1183no rotation is performed when \np{ln\_dynldf\_hor}\forcode{ = .true.}.
1184The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
1185each momentum component.
1186It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.
1187The resulting discrete representation is:
1188\begin{equation}
1189  \label{eq:dyn_ldf_iso}
1190  \begin{split}
1191    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
1192    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_T^{lm}  \left(
1193              {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]
1194                -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}
1195            \right)} \right]}    \right. \\
1196    & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}
1197            }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}
1198            \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}
1199        \right)} \right] \\
1200    &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline
1201              {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }
1202        \right.} \right. \\
1203    &  \ \qquad \qquad \qquad \quad\
1204    - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\
1205    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1206                +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}
1207                \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\
1208    D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\
1209    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_f^{lm}  \left(
1210              {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]
1211                -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}
1212            \right)} \right]}    \right. \\
1213    & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 
1214            }\,\delta_{j} [v] - e_{1t}\, r_{2t}
1215            \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}
1216        \right)} \right] \\
1217    & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 
1218              {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\
1219    &  \ \qquad \qquad \qquad \quad\
1220    - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\
1221    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1222                +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}
1223                \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 
1224  \end{split}
1225\end{equation}
1226where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and
1227the surface of computation ($z$- or $s$-surfaces).
1228The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).
1229
1230%--------------------------------------------------------------------------------------------------------------
1231%           Iso-level bilaplacian operator
1232%--------------------------------------------------------------------------------------------------------------
1233\subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})]
1234            {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})}
1235\label{subsec:DYN_ldf_bilap}
1236
1237The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:dynldf_lap} twice.
1238It requires an additional assumption on boundary conditions:
1239the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
1240while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).
1241%%%
1242\gmcomment{add a remark on the the change in the position of the coefficient}
1243%%%
1244
1245% ================================================================
1246%           Vertical diffusion term
1247% ================================================================
1248\section{Vertical diffusion term (\protect\mdl{dynzdf})}
1249\label{sec:DYN_zdf}
1250%----------------------------------------------namzdf------------------------------------------------------
1251
1252\nlst{namzdf} 
1253%-------------------------------------------------------------------------------------------------------------
1254
1255Options are defined through the \ngn{namzdf} namelist variables.
1256The 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.
1257Two time stepping schemes can be used for the vertical diffusion term:
1258$(a)$ a forward time differencing scheme
1259(\np{ln\_zdfexp}\forcode{ = .true.}) using a time splitting technique (\np{nn\_zdfexp} $>$ 1) or
1260$(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}\forcode{ = .false.})
1261(see \autoref{chap:STP}).
1262Note that namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
1263
1264The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
1265The vertical diffusion operators given by \autoref{eq:PE_zdf} take the following semi-discrete space form:
1266\[
1267  % \label{eq:dynzdf}
1268  \left\{
1269    \begin{aligned}
1270      D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
1271        \ \delta_{k+1/2} [\,u\,]         \right]     \\
1272      \\
1273      D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
1274        \ \delta_{k+1/2} [\,v\,]         \right]
1275    \end{aligned}
1276  \right.
1277\]
1278where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.
1279The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).
1280
1281The surface boundary condition on momentum is the stress exerted by the wind.
1282At the surface, the momentum fluxes are prescribed as the boundary condition on
1283the vertical turbulent momentum fluxes,
1284\begin{equation}
1285  \label{eq:dynzdf_sbc}
1286  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
1287  = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
1288\end{equation}
1289where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in
1290the (\textbf{i},\textbf{j}) coordinate system.
1291The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in
1292the vertical over the mixed layer depth.
1293If the vertical mixing coefficient is small (when no mixed layer scheme is used)
1294the surface stress enters only the top model level, as a body force.
1295The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).
1296
1297The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
1298(see \autoref{sec:ZDF_bfr})
1299
1300% ================================================================
1301% External Forcing
1302% ================================================================
1303\section{External forcings}
1304\label{sec:DYN_forcing}
1305
1306Besides the surface and bottom stresses (see the above section)
1307which are introduced as boundary conditions on the vertical mixing,
1308three other forcings may enter the dynamical equations by affecting the surface pressure gradient.
1309
1310(1) When \np{ln\_apr\_dyn}\forcode{ = .true.} (see \autoref{sec:SBC_apr}),
1311the atmospheric pressure is taken into account when computing the surface pressure gradient.
1312
1313(2) When \np{ln\_tide\_pot}\forcode{ = .true.} and \np{ln\_tide}\forcode{ = .true.} (see \autoref{sec:SBC_tide}),
1314the tidal potential is taken into account when computing the surface pressure gradient.
1315
1316(3) When \np{nn\_ice\_embd}\forcode{ = 2} and LIM or CICE is used
1317($i.e.$ when the sea-ice is embedded in the ocean),
1318the snow-ice mass is taken into account when computing the surface pressure gradient.
1319
1320
1321\gmcomment{ missing : the lateral boundary condition !!!   another external forcing
1322 }
1323
1324% ================================================================
1325% Time evolution term
1326% ================================================================
1327\section{Time evolution term (\protect\mdl{dynnxt})}
1328\label{sec:DYN_nxt}
1329
1330%----------------------------------------------namdom----------------------------------------------------
1331
1332\nlst{namdom} 
1333%-------------------------------------------------------------------------------------------------------------
1334
1335Options are defined through the \ngn{namdom} namelist variables.
1336The general framework for dynamics time stepping is a leap-frog scheme,
1337$i.e.$ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:STP}).
1338The scheme is applied to the velocity, except when
1339using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})
1340in the variable volume case (\key{vvl} defined),
1341where it has to be applied to the thickness weighted velocity (see \autoref{sec:A_momentum}
1342
1343$\bullet$ vector invariant form or linear free surface
1344(\np{ln\_dynhpg\_vec}\forcode{ = .true.} ; \key{vvl} not defined):
1345\[
1346  % \label{eq:dynnxt_vec}
1347  \left\{
1348    \begin{aligned}
1349      &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\
1350      &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
1351    \end{aligned}
1352  \right.
1353\]
1354
1355$\bullet$ flux form and nonlinear free surface
1356(\np{ln\_dynhpg\_vec}\forcode{ = .false.} ; \key{vvl} defined):
1357\[
1358  % \label{eq:dynnxt_flux}
1359  \left\{
1360    \begin{aligned}
1361      &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\
1362      &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
1363      +\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]
1364    \end{aligned}
1365  \right.
1366\]
1367where RHS is the right hand side of the momentum equation,
1368the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.
1369$\gamma$ is initialized as \np{nn\_atfp} (namelist parameter).
1370Its default value is \np{nn\_atfp}\forcode{ = 10.e-3}.
1371In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for
1372the momentum equations.
1373
1374Note that with the filtered free surface,
1375the update of the \textit{after} velocities is done in the \mdl{dynsp\_flt} module,
1376and only array swapping and Asselin filtering is done in \mdl{dynnxt}.
1377
1378% ================================================================
1379\biblio
1380
1381\end{document}
Note: See TracBrowser for help on using the repository browser.