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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Chap_DYN.tex in branches/2017/dev_merge_2017/DOC/TexFiles/Chapters – NEMO

source: branches/2017/dev_merge_2017/DOC/TexFiles/Chapters/Chap_DYN.tex @ 9363

Last change on this file since 9363 was 9363, checked in by nicolasmartin, 6 years ago

Working with htlatex (LaTeX to HTML): disable floatrow package globally, fix several issues by adding \protect to all "fragile" commands (\mdl, \np, \key and \ngn) inside floats (caption or table). Still few specific issues inside Annnex_ISO

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