source: trunk/DOC/TexFiles/Chapters/Chap_DYN.tex @ 2541

Last change on this file since 2541 was 2541, checked in by gm, 10 years ago

v3.3beta: #658 hopefully the final update, including MPP sum, SIGN, IOM, fldread documentation

  • Property svn:executable set to *
File size: 64.8 KB
Line 
1% ================================================================
2% Chapter Ñ Ocean Dynamics (DYN)
3% ================================================================
4\chapter{Ocean Dynamics (DYN)}
5\label{DYN}
6\minitoc
7
8% add a figure for  dynvor ens, ene latices
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 (\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 (\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 \ref{Griffies_al_MWR01, LeclairMadec2009}.
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$.
170Coriolis and momentum advection terms are evaluated using a leapfrog
171scheme, $i.e.$ the velocity appearing in these expressions is centred in
172time (\textit{now} velocity).
173At the lateral boundaries either free slip, no slip or partial slip boundary
174conditions are applied following Chap.\ref{LBC}.
175
176% -------------------------------------------------------------------------------------------------------------
177%        Vorticity term
178% -------------------------------------------------------------------------------------------------------------
179\subsection   [Vorticity term (\textit{dynvor}) ]
180         {Vorticity term (\mdl{dynvor})}
181\label{DYN_vor}
182%------------------------------------------nam_dynvor----------------------------------------------------
183\namdisplay{namdyn_vor} 
184%-------------------------------------------------------------------------------------------------------------
185
186Four discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=true) are available:
187conserving potential enstrophy of horizontally non-divergent flow (ENS scheme) ;
188conserving horizontal kinetic energy (ENE scheme) ; conserving potential enstrophy for
189the relative vorticity term and horizontal kinetic energy for the planetary vorticity
190term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent
191flow and horizontal kinetic energy (ENE scheme) (see  Appendix~\ref{Apdx_C_vor_zad}).
192The vorticity terms are all computed in dedicated routines that can be found in
193the \mdl{dynvor} module.
194
195%-------------------------------------------------------------
196%                 enstrophy conserving scheme
197%-------------------------------------------------------------
198\subsubsection{Enstrophy conserving scheme (\np{ln\_dynvor\_ens}=true)}
199\label{DYN_vor_ens}
200
201In the enstrophy conserving case (ENS scheme), the discrete formulation of the
202vorticity term provides a global conservation of the enstrophy
203($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent
204flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by:
205\begin{equation} \label{Eq_dynvor_ens}
206\left\{ 
207\begin{aligned}
208{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} 
209                                & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2}    \\
210{- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j} 
211                                & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j} 
212\end{aligned} 
213 \right.
214\end{equation} 
215
216%-------------------------------------------------------------
217%                 energy conserving scheme
218%-------------------------------------------------------------
219\subsubsection{Energy conserving scheme (\np{ln\_dynvor\_ene}=true)}
220\label{DYN_vor_ene}
221
222The kinetic energy conserving scheme (ENE scheme) conserves the global
223kinetic energy but not the global enstrophy. It is given by:
224\begin{equation} \label{Eq_dynvor_ene}
225\left\{   \begin{aligned}
226{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
227                            \;  \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
228{- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
229                            \;  \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} }
230\end{aligned}    \right.
231\end{equation} 
232
233%-------------------------------------------------------------
234%                 mix energy/enstrophy conserving scheme
235%-------------------------------------------------------------
236\subsubsection{Mixed energy/enstrophy conserving scheme (\np{ln\_dynvor\_mix}=true) }
237\label{DYN_vor_mix}
238
239For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the
240two previous schemes is used. It consists of the ENS scheme (\ref{Eq_dynvor_ens})
241for the relative vorticity term, and of the ENE scheme (\ref{Eq_dynvor_ene}) applied
242to the planetary vorticity term.
243\begin{equation} \label{Eq_dynvor_mix}
244\left\{ {     \begin{aligned}
245 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} 
246 \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} }
247 \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
248 \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
249{-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j
250 \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} }
251 \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
252 \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill
253\end{aligned}     } \right.
254\end{equation} 
255
256%-------------------------------------------------------------
257%                 energy and enstrophy conserving scheme
258%-------------------------------------------------------------
259\subsubsection{Energy and enstrophy conserving scheme (\np{ln\_dynvor\_een}=true) }
260\label{DYN_vor_een}
261
262In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$ 
263averages of the velocity allows for the presence of grid point oscillation structures
264that will be invisible to the operator. These structures are \textit{computational modes} 
265that will be at least partly damped by the momentum diffusion operator ($i.e.$ the
266subgrid-scale advection), but not by the resolved advection term. The ENS and ENE schemes
267therefore do not contribute to dump any grid point noise in the horizontal velocity field.
268Such noise would result in more noise in the vertical velocity field, an undesirable feature.
269This is a well-known characteristic of $C$-grid discretization where $u$ and $v$ are located
270at different grid points, a price worth paying to avoid a double averaging in the pressure
271gradient term as in the $B$-grid.
272\gmcomment{ To circumvent this, Adcroft (ADD REF HERE)
273Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
274
275A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}.
276The idea is to get rid of the double averaging by considering triad combinations of vorticity.
277It is noteworthy that this solution is conceptually quite similar to the one proposed by
278\citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see App.\ref{Apdx_C}).
279
280The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified
281for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.
282First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
283\begin{equation} \label{Eq_pot_vor}
284q  = \frac{\zeta +f} {e_{3f} }
285\end{equation}
286where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter
287is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:
288\begin{equation} \label{Eq_een_e3f}
289e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
290\end{equation}
291
292%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
293\begin{figure}[!ht]    \begin{center}
294\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_DYN_een_triad.pdf}
295\caption{ \label{Fig_DYN_een_triad} 
296Triads used in the energy and enstrophy conserving scheme (een) for
297$u$-component (upper panel) and $v$-component (lower panel).}
298\end{center}   \end{figure}
299%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
300
301Note that a key point in \eqref{Eq_een_e3f} is that the averaging in the \textbf{i}- and
302\textbf{j}- directions uses the masked vertical scale factor but is always divided by
303$4$, not by the sum of the masks at the four $T$-points. This preserves the continuity of
304$e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and
305extends by continuity the value of $e_{3f}$ into the land areas. This feature is essential for
306the $z$-coordinate with partial steps.
307
308Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
309the following triad combinations of the neighbouring potential vorticities defined at f-points
310(Fig.~\ref{Fig_DYN_een_triad}):
311\begin{equation} \label{Q_triads}
312_i^j \mathbb{Q}^{i_p}_{j_p}
313= \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)
314\end{equation}
315where 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$.
316
317Finally, the vorticity terms are represented as:
318\begin{equation} \label{Eq_dynvor_een}
319\left\{ {
320\begin{aligned}
321 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
322                         {^{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}   \\
323 - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}} 
324                         {^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}   \\
325\end{aligned} 
326} \right.
327\end{equation} 
328
329This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes.
330It conserves both total energy and potential enstrophy in the limit of horizontally
331nondivergent flow ($i.e.$ $\chi$=$0$) (see  Appendix~\ref{Apdx_C_vor_zad}).
332Applied to a realistic ocean configuration, it has been shown that it leads to a significant
333reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.
334Furthermore, used in combination with a partial steps representation of bottom topography,
335it improves the interaction between current and topography, leading to a larger
336topostrophy of the flow  \citep{Barnier_al_OD06, Penduff_al_OS07}.
337
338%--------------------------------------------------------------------------------------------------------------
339%           Kinetic Energy Gradient term
340%--------------------------------------------------------------------------------------------------------------
341\subsection   [Kinetic Energy Gradient term (\textit{dynkeg})]
342         {Kinetic Energy Gradient term (\mdl{dynkeg})}
343\label{DYN_keg}
344
345As demonstrated in Appendix~\ref{Apdx_C}, there is a single discrete formulation
346of the kinetic energy gradient term that, together with the formulation chosen for
347the vertical advection (see below), conserves the total kinetic energy:
348\begin{equation} \label{Eq_dynkeg}
349\left\{ \begin{aligned}
350 -\frac{1}{2 \; e_{1u} }  & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\
351 -\frac{1}{2 \; e_{2v} }  & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   
352\end{aligned} \right.
353\end{equation} 
354
355%--------------------------------------------------------------------------------------------------------------
356%           Vertical advection term
357%--------------------------------------------------------------------------------------------------------------
358\subsection   [Vertical advection term (\textit{dynzad}) ]
359         {Vertical advection term (\mdl{dynzad}) }
360\label{DYN_zad}
361
362The discrete formulation of the vertical advection, together with the formulation
363chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic
364energy. Indeed, the change of KE due to the vertical advection is exactly
365balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}).
366\begin{equation} \label{Eq_dynzad}
367\left\{     \begin{aligned}
368-\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}  \\
369-\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} 
370\end{aligned}         \right.
371\end{equation} 
372
373% ================================================================
374% Coriolis and Advection : flux form
375% ================================================================
376\section{Coriolis and Advection: flux form}
377\label{DYN_adv_cor_flux}
378%------------------------------------------nam_dynadv----------------------------------------------------
379\namdisplay{namdyn_adv} 
380%-------------------------------------------------------------------------------------------------------------
381
382In the flux form (as in the vector invariant form), the Coriolis and momentum
383advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity
384appearing in their expressions is centred in time (\textit{now} velocity). At the
385lateral boundaries either free slip, no slip or partial slip boundary conditions
386are applied following Chap.\ref{LBC}.
387
388
389%--------------------------------------------------------------------------------------------------------------
390%           Coriolis plus curvature metric terms
391%--------------------------------------------------------------------------------------------------------------
392\subsection   [Coriolis plus curvature metric terms (\textit{dynvor}) ]
393         {Coriolis plus curvature metric terms (\mdl{dynvor}) }
394\label{DYN_cor_flux}
395
396In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis
397parameter has been modified to account for the "metric" term. This altered
398Coriolis parameter is thus discretised at $f$-points. It is given by:
399\begin{multline} \label{Eq_dyncor_metric}
400f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right\\
401   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right] 
402                                                                 -  \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right]  }  \ \right)
403\end{multline} 
404
405Any of the (\ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een})
406schemes can be used to compute the product of the Coriolis parameter and the
407vorticity. However, the energy-conserving scheme  (\ref{Eq_dynvor_een}) has
408exclusively been used to date. This term is evaluated using a leapfrog scheme,
409$i.e.$ the velocity is centred in time (\textit{now} velocity).
410
411%--------------------------------------------------------------------------------------------------------------
412%           Flux form Advection term
413%--------------------------------------------------------------------------------------------------------------
414\subsection   [Flux form Advection term (\textit{dynadv}) ]
415         {Flux form Advection term (\mdl{dynadv}) }
416\label{DYN_adv_flux}
417
418The discrete expression of the advection term is given by :
419\begin{equation} \label{Eq_dynadv}
420\left\{ 
421\begin{aligned}
422\frac{1}{e_{1u}\,e_{2u}\,e_{3u}} 
423\left(      \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i       }  \ u_t      \right]   
424          + \delta _{j       } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2}  \ u_f      \right] \right\ \;   \\
425\left.   + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2}  \ u_{uw} \right] \right)   \\
426\\
427\frac{1}{e_{1v}\,e_{2v}\,e_{3v}} 
428\left(     \delta _{i       } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f       \right] 
429         + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i       } \ v_t       \right] \right\ \, \, \\
430\left.  + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw}  \right] \right) \\
431\end{aligned}
432\right.
433\end{equation}
434
435Two advection schemes are available: a $2^{nd}$ order centered finite
436difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS.
437The latter is described in \citet{Shchepetkin_McWilliams_OM05}. The schemes are
438selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.
439In flux form, the schemes differ by the choice of a space and time interpolation to
440define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells,
441$i.e.$ at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and
442$vw$-points for $v$.
443
444%-------------------------------------------------------------
445%                 2nd order centred scheme
446%-------------------------------------------------------------
447\subsubsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_dynadv\_cen2}=true)}
448\label{DYN_adv_cen2}
449
450In the centered $2^{nd}$ order formulation, the velocity is evaluated as the
451mean of the two neighbouring points :
452\begin{equation} \label{Eq_dynadv_cen2}
453\left\{     \begin{aligned}
454 u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\
455 v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j      \quad &  v_{vw}^{cen2} &=\overline v ^{k+1/2}  \\
456\end{aligned}      \right.
457\end{equation} 
458
459The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive
460($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be
461used in conjunction with an explicit diffusion operator to produce a sensible solution.
462The associated time-stepping is performed using a leapfrog scheme in conjunction
463with an Asselin time-filter, so $u$ and $v$ are the \emph{now} velocities.
464
465%-------------------------------------------------------------
466%                 UBS scheme
467%-------------------------------------------------------------
468\subsubsection{Upstream Biased Scheme (UBS) (\np{ln\_dynadv\_ubs}=true)}
469\label{DYN_adv_ubs}
470
471The UBS advection scheme is an upstream biased third order scheme based on
472an upstream-biased parabolic interpolation. For example, the evaluation of
473$u_T^{ubs} $ is done as follows:
474\begin{equation} \label{Eq_dynadv_ubs}
475u_T^{ubs} =\overline u ^i-\;\frac{1}{6}   \begin{cases}
476      u"_{i-1/2}&    \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  \geqslant 0$ }    \\
477      u"_{i+1/2}&    \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  < 0$ }
478\end{cases}
479\end{equation}
480where $u"_{i+1/2} =\delta _{i+1/2} \left[ {\delta _i \left[ u \right]} \right]$. This results
481in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}.
482The overall performance of the advection scheme is similar to that reported in
483\citet{Farrow1995}. It is a relatively good compromise between accuracy and
484smoothness. It is not a \emph{positive} scheme, meaning that false extrema are
485permitted. But the amplitudes of the false extrema are significantly reduced over
486those in the centred second order method. As the scheme already includes
487a diffusion component, it can be used without explicit  lateral diffusion on momentum
488($i.e.$ \np{ln\_dynldf\_lap}=\np{ln\_dynldf\_bilap}=false), and it is recommended to do so.
489
490The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ 
491order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and
492$u_{vw}^{ubs}$ in \eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is
493associated with vertical mixing of momentum. \gmcomment{ gm  pursue the
494sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }
495
496For stability reasons,  the first term in (\ref{Eq_dynadv_ubs}), which corresponds
497to a second order centred scheme, is evaluated using the \textit{now} velocity
498(centred in time), while the second term, which is the diffusion part of the scheme,
499is evaluated using the \textit{before} velocity (forward in time). This is discussed
500by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
501
502Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics)
503schemes only differ by one coefficient. Replacing $1/6$ by $1/8$ in
504(\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
505This option is not available through a namelist parameter, since the $1/6$ coefficient
506is hard coded. Nevertheless it is quite easy to make the substitution in the
507\mdl{dynadv\_ubs} module and obtain a QUICK scheme.
508
509Note also that in the current version of \mdl{dynadv\_ubs}, there is also the
510possibility of using a $4^{th}$ order evaluation of the advective velocity as in
511ROMS. This is an error and should be suppressed soon.
512%%%
513\gmcomment{action :  this have to be done}
514%%%
515
516% ================================================================
517%           Hydrostatic pressure gradient term
518% ================================================================
519\section  [Hydrostatic pressure gradient (\textit{dynhpg})]
520      {Hydrostatic pressure gradient (\mdl{dynhpg})}
521\label{DYN_hpg}
522%------------------------------------------nam_dynhpg---------------------------------------------------
523\namdisplay{namdyn_hpg} 
524%-------------------------------------------------------------------------------------------------------------
525
526The key distinction between the different algorithms used for the hydrostatic
527pressure gradient is the vertical coordinate used, since HPG is a \emph{horizontal} 
528pressure gradient, $i.e.$ computed along geopotential surfaces. As a result, any
529tilt of the surface of the computational levels will require a specific treatment to
530compute the hydrostatic pressure gradient.
531
532The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme,
533$i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$), or
534a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip
535boundary conditions are applied.
536
537%--------------------------------------------------------------------------------------------------------------
538%           z-coordinate with full step
539%--------------------------------------------------------------------------------------------------------------
540\subsection   [$z$-coordinate with full step (\np{ln\_dynhpg\_zco}) ]
541         {$z$-coordinate with full step (\np{ln\_dynhpg\_zco}=true)}
542\label{DYN_hpg_zco}
543
544The hydrostatic pressure can be obtained by integrating the hydrostatic equation
545vertically from the surface. However, the pressure is large at great depth while its
546horizontal gradient is several orders of magnitude smaller. This may lead to large
547truncation errors in the pressure gradient terms. Thus, the two horizontal components
548of the hydrostatic pressure gradient are computed directly as follows:
549
550for $k=km$ (surface layer, $jk=1$ in the code)
551\begin{equation} \label{Eq_dynhpg_zco_surf}
552\left\{ \begin{aligned}
553               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k=km} 
554&= \frac{1}{2} g \   \left. \delta _{i+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
555                  \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k=km} 
556&= \frac{1}{2} g \   \left. \delta _{j+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
557\end{aligned} \right.
558\end{equation} 
559
560for $1<k<km$ (interior layer)
561\begin{equation} \label{Eq_dynhpg_zco}
562\left\{ \begin{aligned}
563               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k} 
564&=             \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k-1} 
565+    \frac{1}{2}\;g\;   \left. \delta _{i+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
566                  \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k} 
567&=                \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k-1} 
568+    \frac{1}{2}\;g\;   \left. \delta _{j+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
569\end{aligned} \right.
570\end{equation} 
571
572Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of
573the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface
574level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the
575surface pressure gradient is included in \eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco} 
576through the space and time variations of the vertical scale factor $e_{3w}$.
577
578%--------------------------------------------------------------------------------------------------------------
579%           z-coordinate with partial step
580%--------------------------------------------------------------------------------------------------------------
581\subsection   [$z$-coordinate with partial step (\np{ln\_dynhpg\_zps})]
582         {$z$-coordinate with partial step (\np{ln\_dynhpg\_zps}=true)}
583\label{DYN_hpg_zps}
584
585With partial bottom cells, tracers in horizontally adjacent cells generally live at
586different depths. Before taking horizontal gradients between these tracer points,
587a linear interpolation is used to approximate the deeper tracer as if it actually lived
588at the depth of the shallower tracer point.
589
590Apart from this modification, the horizontal hydrostatic pressure gradient evaluated
591in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case.
592As explained in detail in section \S\ref{TRA_zpshde}, the nonlinearity of pressure
593effects in the equation of state is such that it is better to interpolate temperature and
594salinity vertically before computing the density. Horizontal gradients of temperature
595and salinity are needed for the TRA modules, which is the reason why the horizontal
596gradients of density at the deepest model level are computed in module \mdl{zpsdhe} 
597located in the TRA directory and described in \S\ref{TRA_zpshde}.
598
599%--------------------------------------------------------------------------------------------------------------
600%           s- and s-z-coordinates
601%--------------------------------------------------------------------------------------------------------------
602\subsection{$s$- and $z$-$s$-coordinates}
603\label{DYN_hpg_sco}
604
605Pressure gradient formulations in an $s$-coordinate have been the subject of a vast
606number of papers ($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}).
607A number of different pressure gradient options are coded, but they are not yet fully
608documented or tested.
609
610$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true,
611\np{ln\_dynhpg\_hel}=true)
612\begin{equation} \label{Eq_dynhpg_sco}
613\left\{ \begin{aligned}
614 - \frac{1}                   {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  p^h  \right] 
615+ \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  z_t   \right]    \\
616 - \frac{1}                   {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  p^h  \right] 
617+ \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  z_t   \right]    \\
618\end{aligned} \right.
619\end{equation} 
620
621Where the first term is the pressure gradient along coordinates, computed as in
622\eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of
623the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
624($e_{3w}$). The version \np{ln\_dynhpg\_hel}=true has been added by Aike
625Beckmann and involves a redefinition of the relative position of $T$-points relative
626to $w$-points.
627
628$\bullet$ Weighted density Jacobian (WDJ) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=true)
629
630$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} 
631(\np{ln\_dynhpg\_djc}=true)
632
633$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true)
634
635Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume
636formulation is activated (\key{vvl}) because in that case, even with a flat bottom,
637the coordinate surfaces are not horizontal but follow the free surface
638\citep{Levier2007}. The other pressure gradient options are not yet available.
639
640%--------------------------------------------------------------------------------------------------------------
641%           Time-scheme
642%--------------------------------------------------------------------------------------------------------------
643\subsection   [Time-scheme (\np{ln\_dynhpg\_imp}) ]
644         {Time-scheme (\np{ln\_dynhpg\_imp}= true/false)}
645\label{DYN_hpg_imp}
646
647The default time differencing scheme used for the horizontal pressure gradient is
648a leapfrog scheme and therefore the density used in all discrete expressions given
649above is the  \textit{now} density, computed from the \textit{now} temperature and
650salinity. In some specific cases (usually high resolution simulations over an ocean
651domain which includes weakly stratified regions) the physical phenomenon that
652controls the time-step is internal gravity waves (IGWs). A semi-implicit scheme for
653doubling the stability limit associated with IGWs can be used \citep{Brown_Campana_MWR78,
654Maltrud1998}. It involves the evaluation of the hydrostatic pressure gradient as an
655average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ ($i.e.$ 
656\textit{before}\textit{now} and  \textit{after} time-steps), rather than at the central
657time level $t$ only, as in the standard leapfrog scheme.
658
659$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}=true):
660
661\begin{equation} \label{Eq_dynhpg_lf}
662\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
663   -\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right]
664\end{equation}
665
666$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}=true):
667\begin{equation} \label{Eq_dynhpg_imp}
668\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
669   -\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]
670\end{equation}
671
672The semi-implicit time scheme \eqref{Eq_dynhpg_imp} is made possible without
673significant additional computation since the density can be updated to time level
674$t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can
675be easily shown that the stability limit associated with the hydrostatic pressure
676gradient doubles using \eqref{Eq_dynhpg_imp} compared to that using the
677standard leapfrog scheme \eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp} 
678is equivalent to applying a time filter to the pressure gradient to eliminate high
679frequency IGWs. Obviously, when using \eqref{Eq_dynhpg_imp}, the doubling of
680the time-step is achievable only if no other factors control the time-step, such as
681the stability limits associated with advection or diffusion.
682
683In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=true.
684In this case, we choose to apply the time filter to temperature and salinity used in
685the equation of state, instead of applying it to the hydrostatic pressure or to the
686density, so that no additional storage array has to be defined. The density used to
687compute the hydrostatic pressure gradient (whatever the formulation) is evaluated
688as follows:
689\begin{equation} \label{Eq_rho_flt}
690   \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
691 \quad     \text{with}  \quad 
692   \widetilde{X} = 1 / 4 \left(  X^{t+\rdt} +2 \,X^t + X^{t-\rdt}  \right)
693\end{equation}
694
695Note that in the semi-implicit case, it is necessary to save the filtered density, an
696extra three-dimensional field, in the restart file to restart the model with exact
697reproducibility. This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter.
698
699% ================================================================
700% Surface Pressure Gradient
701% ================================================================
702\section  [Surface pressure gradient (\textit{dynspg}) ]
703      {Surface pressure gradient (\mdl{dynspg})}
704\label{DYN_spg}
705%-----------------------------------------nam_dynspg----------------------------------------------------
706\namdisplay{namdyn_spg} 
707%------------------------------------------------------------------------------------------------------------
708
709$\ $\newline      %force an empty line
710
711%%%
712The surface pressure gradient term is related to the representation of the free surface (\S\ref{PE_hor_pg}). The main distinction is between the fixed volume case (linear free surface) and the variable volume case (nonlinear free surface, \key{vvl} is defined). In the linear free surface case (\S\ref{PE_free_surface}) the vertical scale factors $e_{3}$ are fixed in time, while they are time-dependent in the nonlinear case (\S\ref{PE_free_surface}). With both linear and nonlinear free surface, external gravity waves are allowed in the equations, which imposes a very small time step when an explicit time stepping is used. Two methods are proposed to allow a longer time step for the three-dimensional equations: the filtered free surface, which is a modification of the continuous equations (see \eqref{Eq_PE_flt}), and the split-explicit free surface described below. The extra term introduced in the filtered method is calculated implicitly, so that the update of the next velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
713
714%%%
715
716
717The form of the surface pressure gradient term depends on how the user wants to handle
718the fast external gravity waves that are a solution of the analytical equation (\S\ref{PE_hor_pg}).
719Three formulations are available, all controlled by a CPP key (ln\_dynspg\_xxx):
720an explicit formulation which requires a small time step ;
721a filtered free surface formulation which allows a larger time step by adding a filtering
722term into the momentum equation ;
723and a split-explicit free surface formulation, described below, which also allows a larger time step.
724
725The extra term introduced in the filtered method is calculated
726implicitly, so that a solver is used to compute it. As a consequence the update of the $next$ 
727velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.
728
729
730
731%--------------------------------------------------------------------------------------------------------------
732% Explicit free surface formulation
733%--------------------------------------------------------------------------------------------------------------
734\subsection{Explicit free surface (\key{dynspg\_exp})}
735\label{DYN_spg_exp}
736
737In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step
738is chosen to be small enough to resolve the external gravity waves (typically a few tens of seconds).
739The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time),
740is thus simply given by :
741\begin{equation} \label{Eq_dynspg_exp}
742\left\{ \begin{aligned}
743 - \frac{1}{e_{1u}\,\rho_o} \;   \delta _{i+1/2} \left[  \,\rho \,\eta\,  \right]   \\
744 - \frac{1}{e_{2v}\,\rho_o} \;   \delta _{j+1/2} \left[  \,\rho \,\eta\,  \right] 
745\end{aligned} \right.
746\end{equation} 
747
748Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined), the surface pressure
749gradient is already included in the momentum tendency  through the level thickness variation
750allowed in the computation of the hydrostatic pressure gradient. Thus, nothing is done in the \mdl{dynspg\_exp} module.
751
752%--------------------------------------------------------------------------------------------------------------
753% Split-explict free surface formulation
754%--------------------------------------------------------------------------------------------------------------
755\subsection{Split-Explicit free surface (\key{dynspg\_ts})}
756\label{DYN_spg_ts}
757
758The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),
759also called the time-splitting formulation, follows the one
760proposed by \citet{Griffies_Bk04}. The general idea is to solve the free surface
761equation and the associated barotropic velocity equations with a smaller time
762step than $\rdt$, the time step used for the three dimensional prognostic
763variables (Fig.~\ref{Fig_DYN_dynspg_ts}).
764The size of the small time step, $\rdt_e$ (the external mode or barotropic time step)
765 is provided through the \np{nn\_baro} namelist parameter as:
766$\rdt_e = \rdt / nn\_baro$.
767 
768
769%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
770\begin{figure}[!t]    \begin{center}
771\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf}
772\caption{  \label{Fig_DYN_dynspg_ts}
773Schematic of the split-explicit time stepping scheme for the external
774and internal modes. Time increases to the right.
775Internal mode time steps (which are also the model time steps) are denoted
776by $t-\rdt$, $t, t+\rdt$, and $t+2\rdt$.
777The curved line represents a leap-frog time step, and the smaller time
778steps $N \rdt_e=\frac{3}{2}\rdt$ are denoted by the zig-zag line.
779The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ 
780represents the interaction between the external and internal motions.
781While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog
782integration carries the external mode variables (surface height and vertically
783integrated velocity) from $t$ to $t+\frac{3}{2} \rdt$ using N external time
784steps of length $\rdt_e$. Time averaging the external fields over the
785$\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated
786velocity and the sea surface height at the model timestep $t+\rdt$.
787These averaged values are used to update \textbf{M}(t) with both the surface
788pressure gradient and the Coriolis force, therefore providing the $t+\rdt$
789velocity.  The model time stepping scheme can then be achieved by a baroclinic
790leap-frog time step that carries the surface height from $t-\rdt$ to $t+\rdt$}
791\end{center}    \end{figure}
792%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
793
794The split-explicit formulation has a damping effect on external gravity waves,
795which is weaker damping than that for the filtered free surface but still significant, as
796shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.
797
798%>>>>>===============
799\gmcomment{               %%% copy from griffies Book
800
801\textbf{title: Time stepping the barotropic system }
802
803Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$. Hence,
804we can update the surface height and vertically integrated velocity with a leap-frog
805scheme using the small barotropic time step $\rdt$. We have
806
807\begin{equation} \label{DYN_spg_ts_eta}
808\eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
809   = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] 
810\end{equation}
811\begin{multline} \label{DYN_spg_ts_u}
812\textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1}\\
813   = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n})
814   - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right]
815\end{multline}
816\
817
818In 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)}$ 
819and $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
820that sets the barotropic time steps via
821\begin{equation} \label{DYN_spg_ts_t}
822t_n=\tau+n\rdt   
823\end{equation}
824with $n$ an integer. The density scaled surface pressure is evaluated via
825\begin{equation} \label{DYN_spg_ts_ps}
826p_s^{(b)}(\tau,t_{n}) = \begin{cases}
827   g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_&      \text{non-linear case} \\
828   g \;\eta_s^{(b)}(\tau,t_{n}&      \text{linear case} 
829   \end{cases}
830\end{equation}
831To get started, we assume the following initial conditions
832\begin{equation} \label{DYN_spg_ts_eta}
833\begin{split}
834\eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)}
835\\
836\eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} 
837\end{split}
838\end{equation}
839with
840\begin{equation} \label{DYN_spg_ts_etaF}
841 \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
842\end{equation}
843the time averaged surface height taken from the previous barotropic cycle. Likewise,
844\begin{equation} \label{DYN_spg_ts_u}
845\textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)}   \\
846\\
847\textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}   
848\end{equation}
849with
850\begin{equation} \label{DYN_spg_ts_u}
851 \overline{\textbf{U}^{(b)}(\tau)} 
852   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n})
853\end{equation}
854the time averaged vertically integrated transport. Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration.
855
856Upon 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$ 
857\begin{equation} \label{DYN_spg_ts_u}
858\textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} 
859   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n})
860\end{equation}
861The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using the following form
862
863\begin{equation} \label{DYN_spg_ts_ssh}
864\eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right] 
865\end{equation}
866
867 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).
868 
869In general, some form of time filter is needed to maintain integrity of the surface
870height field due to the leap-frog splitting mode in equation \ref{DYN_spg_ts_ssh}. We
871have tried various forms of such filtering, with the following method discussed in
872\cite{Griffies_al_MWR01} chosen due to its stability and reasonably good maintenance of
873tracer conservation properties (see Section ??)
874
875\begin{equation} \label{DYN_spg_ts_sshf}
876\eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)} 
877\end{equation}
878Another approach tried was
879
880\begin{equation} \label{DYN_spg_ts_sshf2}
881\eta^{F}(\tau-\Delta) = \eta(\tau)
882   + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
883                + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
884\end{equation}
885
886which is useful since it isolates all the time filtering aspects into the term multiplied
887by $\alpha$. This isolation allows for an easy check that tracer conservation is exact when
888eliminating 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.
889
890}            %%end gm comment (copy of griffies book)
891
892%>>>>>===============
893
894
895%--------------------------------------------------------------------------------------------------------------
896% Filtered free surface formulation
897%--------------------------------------------------------------------------------------------------------------
898\subsection{Filtered free surface (\key{dynspg\_flt})}
899\label{DYN_spg_fltp}
900
901The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.
902The extra term introduced in the equations (see \S\ref{PE_free_surface}) is solved implicitly.
903The elliptic solvers available in the code are documented in \S\ref{MISC}.
904
905%% gm %%======>>>>   given here the discrete eqs provided to the solver
906\gmcomment{               %%% copy from chap-model basics
907\begin{equation} \label{Eq_spg_flt}
908\frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
909- g \nabla \left( \tilde{\rho} \ \eta \right)
910- g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)
911\end{equation}
912where $T_c$, is a parameter with dimensions of time which characterizes the force,
913$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ 
914represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
915non-linear and viscous terms in \eqref{Eq_PE_dyn}.
916}   %end gmcomment
917
918Note that in the linear free surface formulation (\key{vvl} not defined), the ocean depth
919is time-independent and so is the matrix to be inverted. It is computed once and for all and applies to all ocean time steps.
920
921% ================================================================
922% Lateral diffusion term
923% ================================================================
924\section  [Lateral diffusion term (\textit{dynldf})]
925      {Lateral diffusion term (\mdl{dynldf})}
926\label{DYN_ldf}
927%------------------------------------------nam_dynldf----------------------------------------------------
928\namdisplay{namdyn_ldf} 
929%-------------------------------------------------------------------------------------------------------------
930
931The options available for lateral diffusion are to use either laplacian
932(rotated or not) or biharmonic operators. The coefficients may be constant
933or spatially variable; the description of the coefficients is found in the chapter
934on lateral physics (Chap.\ref{LDF}). The lateral diffusion of momentum is
935evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression
936is the \textit{before} velocity in time, except for the pure vertical component
937that appears when a tensor of rotation is used. This latter term is solved
938implicitly together with the vertical diffusion term (see \S\ref{DOM_nxt})
939
940At the lateral boundaries either free slip, no slip or partial slip boundary
941conditions are applied according to the user's choice (see Chap.\ref{LBC}).
942
943% ================================================================
944\subsection   [Iso-level laplacian operator (\np{ln\_dynldf\_lap}) ]
945         {Iso-level laplacian operator (\np{ln\_dynldf\_lap}=true)}
946\label{DYN_ldf_lap}
947
948For lateral iso-level diffusion, the discrete operator is:
949\begin{equation} \label{Eq_dynldf_lap}
950\left\{ \begin{aligned}
951 D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm} 
952\;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta _j \left[
953{A_f^{lm} \;e_{3f} \zeta } \right] \\ 
954\\
955 D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta _{j+1/2} \left[ {A_T^{lm} 
956\;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta _i \left[
957{A_f^{lm} \;e_{3f} \zeta } \right] \\ 
958\end{aligned} \right.
959\end{equation} 
960
961As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence
962and curl of the vorticity) preserves symmetry and ensures a complete
963separation between the vorticity and divergence parts of the momentum diffusion.
964
965%--------------------------------------------------------------------------------------------------------------
966%           Rotated laplacian operator
967%--------------------------------------------------------------------------------------------------------------
968\subsection   [Rotated laplacian operator (\np{ln\_dynldf\_iso}) ]
969         {Rotated laplacian operator (\np{ln\_dynldf\_iso}=true)}
970\label{DYN_ldf_iso}
971
972A rotation of the lateral momentum diffusion operator is needed in several cases:
973for iso-neutral diffusion in the $z$-coordinate (\np{ln\_dynldf\_iso}=true) and for
974either iso-neutral (\np{ln\_dynldf\_iso}=true) or geopotential
975(\np{ln\_dynldf\_hor}=true) diffusion in the $s$-coordinate. In the partial step
976case, coordinates are horizontal except at the deepest level and no
977rotation is performed when \np{ln\_dynldf\_hor}=true. The diffusion operator
978is defined simply as the divergence of down gradient momentum fluxes on each
979momentum component. It must be emphasized that this formulation ignores
980constraints on the stress tensor such as symmetry. The resulting discrete
981representation is:
982\begin{equation} \label{Eq_dyn_ldf_iso}
983\begin{split}
984 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
985&  \left\{\quad  {\delta _{i+1/2} \left[ {A_T^{lm}  \left(
986    {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta _{i}[u]
987   -e_{2t} \; r_{1t} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}}
988 \right)} \right]}   \right.
989\\ 
990& \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} 
991}\,\delta _{j+1/2} [u] - e_{1f}\, r_{2f} 
992\,\overline{\overline {\delta _{k+1/2} [u]}} ^{\,j+1/2,\,k}} 
993\right)} \right]
994\\ 
995&\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline 
996{\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } 
997\right.} \right.
998\\ 
999&  \ \qquad \qquad \qquad \quad\
1000- e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2}
1001\\ 
1002& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1003+\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} 
1004\right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} 
1005\\
1006\\
1007 D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} }    \\
1008&  \left\{\quad  {\delta _{i+1/2} \left[ {A_f^{lm}  \left(
1009    {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta _{i+1/2}[v]
1010   -e_{2f} \; r_{1f} \,\overline{\overline {\delta _{k+1/2}[v]}}^{\,i+1/2,\,k}}
1011 \right)} \right]}   \right.
1012\\ 
1013& \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 
1014}\,\delta _{j} [v] - e_{1t}\, r_{2t} 
1015\,\overline{\overline {\delta _{k+1/2} [v]}} ^{\,j,\,k}} 
1016\right)} \right]
1017\\ 
1018& \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 
1019{\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right.
1020\\
1021&  \ \qquad \qquad \qquad \quad\
1022- e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2}
1023\\ 
1024& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
1025+\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} 
1026\right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 
1027 \end{split}
1028\end{equation}
1029where $r_1$ and $r_2$ are the slopes between the surface along which the
1030diffusion operator acts and the surface of computation ($z$- or $s$-surfaces).
1031The way these slopes are evaluated is given in the lateral physics chapter
1032(Chap.\ref{LDF}).
1033
1034%--------------------------------------------------------------------------------------------------------------
1035%           Iso-level bilaplacian operator
1036%--------------------------------------------------------------------------------------------------------------
1037\subsection   [Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap})]
1038         {Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap}=true)}
1039\label{DYN_ldf_bilap}
1040
1041The lateral fourth order operator formulation on momentum is obtained by
1042applying \eqref{Eq_dynldf_lap} twice. It requires an additional assumption on
1043boundary conditions: the first derivative term normal to the coast depends on
1044the free or no-slip lateral boundary conditions chosen, while the third
1045derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}).
1046%%%
1047\gmcomment{add a remark on the the change in the position of the coefficient}
1048%%%
1049
1050% ================================================================
1051%           Vertical diffusion term
1052% ================================================================
1053\section  [Vertical diffusion term (\mdl{dynzdf})]
1054      {Vertical diffusion term (\mdl{dynzdf})}
1055\label{DYN_zdf}
1056%----------------------------------------------namzdf------------------------------------------------------
1057\namdisplay{namzdf} 
1058%-------------------------------------------------------------------------------------------------------------
1059
1060The large vertical diffusion coefficient found in the surface mixed layer together
1061with high vertical resolution implies that in the case of explicit time stepping there
1062would be too restrictive a constraint on the time step. Two time stepping schemes
1063can be used for the vertical diffusion term : $(a)$ a forward time differencing
1064scheme (\np{ln\_zdfexp}=true) using a time splitting technique
1065(\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme
1066(\np{ln\_zdfexp}=false) (see \S\ref{DOM_nxt}). Note that namelist variables
1067\np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
1068
1069The formulation of the vertical subgrid scale physics is the same whatever
1070the vertical coordinate is. The vertical diffusion operators given by
1071\eqref{Eq_PE_zdf} take the following semi-discrete space form:
1072\begin{equation} \label{Eq_dynzdf}
1073\left\{   \begin{aligned}
1074D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
1075                              \ \delta _{k+1/2} [\,u\,]         \right]     \\
1076\\
1077D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta _k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
1078                              \ \delta _{k+1/2} [\,v\,]         \right]
1079\end{aligned}   \right.
1080\end{equation} 
1081where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and
1082diffusivity coefficients. The way these coefficients are evaluated
1083depends on the vertical physics used (see \S\ref{ZDF}).
1084
1085The surface boundary condition on momentum is the stress exerted by
1086the wind. At the surface, the momentum fluxes are prescribed as the boundary
1087condition on the vertical turbulent momentum fluxes,
1088\begin{equation} \label{Eq_dynzdf_sbc}
1089\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
1090    = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v }
1091\end{equation}
1092where $\left( \tau _u ,\tau _v \right)$ are the two components of the wind stress
1093vector in the (\textbf{i},\textbf{j}) coordinate system. The high mixing coefficients
1094in the surface mixed layer ensure that the surface wind stress is distributed in
1095the vertical over the mixed layer depth. If the vertical mixing coefficient
1096is small (when no mixed layer scheme is used) the surface stress enters only
1097the top model level, as a body force. The surface wind stress is calculated
1098in the surface module routines (SBC, see Chap.\ref{SBC})
1099
1100The turbulent flux of momentum at the bottom of the ocean is specified through
1101a bottom friction parameterisation (see \S\ref{ZDF_bfr})
1102
1103% ================================================================
1104% External Forcing
1105% ================================================================
1106\section{External Forcings}
1107\label{DYN_forcing}
1108
1109Besides the surface and bottom stresses (see the above section) which are
1110introduced as boundary conditions on the vertical mixing, two other forcings
1111enter the dynamical equations.
1112
1113One is the effect of atmospheric pressure on the ocean dynamics.
1114Another forcing term is the tidal potential.
1115Both of which will be introduced into the reference version soon.
1116
1117\gmcomment{atmospheric pressure is there!!!!    include its description }
1118
1119% ================================================================
1120% Time evolution term
1121% ================================================================
1122\section  [Time evolution term (\textit{dynnxt})]
1123      {Time evolution term (\mdl{dynnxt})}
1124\label{DYN_nxt}
1125
1126%----------------------------------------------namdom----------------------------------------------------
1127\namdisplay{namdom} 
1128%-------------------------------------------------------------------------------------------------------------
1129
1130The general framework for dynamics time stepping is a leap-frog scheme,
1131$i.e.$ a three level centred time scheme associated with an Asselin time filter
1132(cf. Chap.\ref{STP}). The scheme is applied to the velocity, except when using
1133the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in the variable
1134volume case (\key{vvl} defined), where it has to be applied to the thickness
1135weighted velocity (see \S\ref{Apdx_A_momentum}
1136
1137$\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}=true ; \key{vvl} not defined):
1138\begin{equation} \label{Eq_dynnxt_vec}
1139\left\{   \begin{aligned}
1140&u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\
1141&u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
1142\end{aligned}   \right.
1143\end{equation} 
1144
1145$\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}=false ; \key{vvl} defined):
1146\begin{equation} \label{Eq_dynnxt_flux}
1147\left\{   \begin{aligned}
1148&\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\
1149&\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
1150  +\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]
1151\end{aligned}   \right.
1152\end{equation} 
1153where RHS is the right hand side of the momentum equation, the subscript $f$ 
1154denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is
1155initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = $10^{-3}$.
1156In both cases, the modified Asselin filter is not applied since perfect conservation
1157is not an issue for the momentum equations.
1158
1159Note that with the filtered free surface, the update of the \textit{after} velocities
1160is done in the \mdl{dynsp\_flt} module, and only array swapping
1161and Asselin filtering is done in \mdl{dynnxt}.
1162
1163% ================================================================
Note: See TracBrowser for help on using the repository browser.