# 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/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters – NEMO

# source:branches/2015/nemo_v3_6_STABLE/DOC/TexFiles/Chapters/Chap_DYN.tex@6317

Last change on this file since 6317 was 6317, checked in by mathiot, 8 years ago

ISF: update documentation and biblio

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