source: branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex @ 2197

% ================================================================
% Chapter Ñ Ocean Dynamics (DYN)
% ================================================================
\chapter{Ocean Dynamics (DYN)}
\label{DYN}
\minitoc

% add a figure for  dynvor ens, ene latices

%\vspace{2.cm}
$\ $\newline      %force an empty line

Using the representation described in Chap.\ref{DOM}, several semi-discrete 
space forms of the dynamical equations are available depending on the vertical 
coordinate used and on the conservation properties of the vorticity term. In all 
the equations presented here, the masking has been omitted for simplicity. 
One must be aware that all the quantities are masked fields and that each time a 
average or difference operator is used, the resulting field is multiplied by a mask.

The prognostic ocean dynamics equation can be summarized as follows:
\begin{equation*}
\text{NXT} = \dbinom {\text{VOR} + \text{KEG} + \text {ZAD} }
                  {\text{COR} + \text{ADV}                       }
         + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
\end{equation*}
NXT stands for next, referring to the time-stepping. The first group of terms on 
the rhs of the momentum equations corresponds to the Coriolis and advection 
terms that are decomposed into a vorticity part (VOR), a kinetic energy part (KEG) 
and, a vertical advection part (ZAD) in the vector invariant formulation or a Coriolis 
and advection part (COR+ADV) in the flux formulation. The terms following these 
are the pressure gradient contributions (HPG, Hydrostatic Pressure Gradient, 
and SPG, Surface Pressure Gradient); and contributions from lateral diffusion 
(LDF) and vertical diffusion (ZDF), which are added to the rhs in the \mdl{dynldf} 
and \mdl{dynzdf} modules. The vertical diffusion term includes the surface and 
bottom stresses. The external forcings and parameterisations require complex 
inputs (surface wind stress calculation using bulk formulae, estimation of mixing 
coefficients) that are carried out in modules SBC, LDF and ZDF and are described 
in Chapters \ref{SBC}, \ref{LDF} and \ref{ZDF}, respectively. 

In the present chapter we also describe the diagnostic equations used to compute 
the horizontal divergence and curl of the velocities (\emph{divcur} module) as well 
as the vertical velocity (\emph{wzvmod} module).

The different options available to the user are managed by namelist variables. 
For equation term \textit{ttt}, the logical namelist variables are \textit{ln\_dynttt\_xxx}, 
where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 
If a CPP key is used for this term its name is \textbf{key\_ttt}. The corresponding 
code can be found in the \textit{dynttt\_xxx} module in the DYN directory, and it is 
usually computed in the \textit{dyn\_ttt\_xxx} subroutine.

The user has the option of extracting each tendency term of both the rhs of the 
3D momentum equation (\key{trddyn} defined) for output, as described in 
Chap.\ref{MISC}.  Furthermore, the tendency terms associated to the 2D 
barotropic vorticity balance (\key{trdvor} defined) can be derived on-line from the 
3D terms.
%%%
\gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does 
MISC correspond to "extracting tendency terms" or "vorticity balance"?}

$\ $\newline    % force a new ligne

% ================================================================
% Sea Surface Height evolution & Diagnostics variables
% ================================================================
\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)}
\label{DYN_divcur_wzv}

%--------------------------------------------------------------------------------------------------------------
%           Horizontal divergence and relative vorticity
%--------------------------------------------------------------------------------------------------------------
\subsection   [Horizontal divergence and relative vorticity (\textit{divcur})]
         {Horizontal divergence and relative vorticity (\mdl{divcur})}
\label{DYN_divcur}

The vorticity is defined at an $f$-point ($i.e.$ corner point) as follows:
\begin{equation} \label{Eq_divcur_cur}
\zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta _{i+1/2} \left[ {e_{2v}\;v} \right]
                          -\delta _{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right)
\end{equation} 

The horizontal divergence is defined at a $T$-point. It is given by:
\begin{equation} \label{Eq_divcur_div}
\chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
      \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u} \right]
             +\delta _j \left[ {e_{1v}\,e_{3v}\,v} \right]} \right)
\end{equation} 

Note that in the $z$-coordinate with full step (\key{zco} is defined), 
$e_{3u}$=$e_{3v}$=$e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}.

Note also that whereas the vorticity have the same discrete expression in $z$- 
and $s$-coordinate, its physical meaning is not identical. $\zeta$ is a pseudo 
vorticity along $s$-surfaces (only pseudo because $(u,v)$ are still defined along 
geopotential surfaces, but are no more necessary defined at the same depth).

The vorticity and divergence at the \textit{before} step are used in the computation 
of the horizontal diffusion of momentum. Note that because they have been 
calculated prior to the Asselin filtering of the \textit{before} velocities, the 
\textit{before} vorticity and divergence arrays must be included in the restart file 
to ensure perfect restartability. The vorticity and divergence at the \textit{now} 
time step are used for the computation of the nonlinear advection and of the 
vertical velocity respectively. 

%--------------------------------------------------------------------------------------------------------------
%           Sea Surface Height evolution
%--------------------------------------------------------------------------------------------------------------
\subsection   [Sea surface height evolution and vertical velocity (\textit{sshwzv})]
         {Horizontal divergence and relative vorticity (\mdl{sshwzv})}
\label{DYN_sshwzv}

The sea surface height is given by :
\begin{equation} \label{Eq_dynspg_ssh}
\begin{aligned}
\frac{\partial \eta }{\partial t}
&\equiv    \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left(  \delta _i \left[ {e_{2u}\,e_{3u}\;u} \right]
                                                                                  +\delta _j \left[ {e_{1v}\,e_{3v}\;v} \right]  \right) } 
           -    \frac{\textit{emp}}{\rho _w }   \\
&\equiv    \sum\limits_k {\chi \ e_{3t}}  -  \frac{\textit{emp}}{\rho _w }
\end{aligned}
\end{equation}
where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), 
expressed in Kg/m$^2$/s (which is equal to mm/s), and $\rho _w$=1,000~Kg/m$^3$ 
is the volumic mass of pure water. If river runoff is expressed as a surface freshwater 
flux (see \S\ref{SBC}) then \textit{emp} can be written as the evaporation minus 
precipitation, minus the river runoff. The sea-surface height is evaluated 
using exactly the same time stepping as the tracer equation \eqref{Eq_tra_nxt}: 
a leapfrog scheme in combination with an Asselin time filter, $i.e.$ the velocity appearing 
in \eqref{Eq_dynspg_ssh} is centred in time (\textit{now} velocity). 
This is of paramount importance. Substituing $T$ by $1$ in the tracer equation and summing
over the water column must lead to the sea surface height equation otherwise tracer content
could not be conserved \ref{Griffies_al_MWR01, LeclairMadec2009}.

The vertical velocity is computed by an upward integration of the horizontal 
divergence from the bottom, taken into account the change of the thickness of the levels :

\begin{equation} \label{Eq_wzv}
\left\{   \begin{aligned}
&\left. w \right|_{3/2} \quad= 0    \\
&\left. w \right|_{k+1/2}     = \left. w \right|_{k+1/2}  + e_{3t}\;  \left. \chi \right|_k  
                                         - \frac{ e_{3t}^{t+1} - e_{3t}^{t-1} } {2 \rdt}
\end{aligned}   \right.
\end{equation}

In case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$, 
as changes in the the divergence of the barotropic transport is absorbed in the change 
of the levels thickness.re-oriented downward co
In case of linear free surface, the time derivative in \eqref{Eq_wzv} cancel out.
The upper boundary condition applies at a fixed level $z=0$. The top vertical velocity 
is thus equal to the divergence of the barotropic transport ($i.e.$ the first term in the
right-hand-side of \eqref{Eq_dynspg_ssh}).

Note also that whereas the vertical velocity has the same discrete 
expression in $z$- and $s$-coordinate, its physical meaning is not the same: 
in the second case, $w$ is the velocity normal to the $s$-surfaces. 
Note also that the $k$-axis is re-oriented downward in the \textsc{fortran} code compare 
to the indexing used in the semi-discrete equations such as \eqref{Eq_wzv} 
(see  \S\ref{DOM_Num_Index_vertical}). 


% ================================================================
% Coriolis and Advection terms: vector invariant form
% ================================================================
\section{Coriolis and Advection: vector invariant form}
\label{DYN_adv_cor_vect}
%-----------------------------------------nam_dynadv----------------------------------------------------
\namdisplay{namdyn_adv} 
%-------------------------------------------------------------------------------------------------------------

The vector invariant form of the momentum equations is the one most 
often used in applications of \NEMO ocean model. The flux form option 
(see next section) has been introduced since version $2$. 
Coriolis and momentum advection terms are evaluated using a leapfrog 
scheme, $i.e.$ the velocity appearing in these expressions is centred in 
time (\textit{now} velocity). 
At the lateral boundaries either free slip, no slip or partial slip boundary 
conditions are applied following Chap.\ref{LBC}.

% -------------------------------------------------------------------------------------------------------------
%        Vorticity term 
% -------------------------------------------------------------------------------------------------------------
\subsection   [Vorticity term (\textit{dynvor}) ]
         {Vorticity term (\mdl{dynvor})}
\label{DYN_vor}
%------------------------------------------nam_dynvor----------------------------------------------------
\namdisplay{namdyn_vor} 
%-------------------------------------------------------------------------------------------------------------

Four discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=true) are available: 
conserving potential enstrophy of horizontally non-divergent flow (ENS scheme) ; 
conserving horizontal kinetic energy (ENE scheme) ; conserving potential enstrophy for 
the relative vorticity term and horizontal kinetic energy for the planetary vorticity 
term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent 
flow and horizontal kinetic energy (ENE scheme) (see  Appendix~\ref{Apdx_C_vor_zad}). 
The vorticity terms are given below for the general case, but note that in the full step 
$z$-coordinate (\key{zco} is defined), $e_{3u}$=$e_{3v}$=$e_{3f}$ so that the vertical scale 
factors disappear. They are all computed in dedicated routines that can be found in 
the \mdl{dynvor} module.

%-------------------------------------------------------------
%                 enstrophy conserving scheme
%-------------------------------------------------------------
\subsubsection{Enstrophy conserving scheme (\np{ln\_dynvor\_ens}=true)}
\label{DYN_vor_ens}

In the enstrophy conserving case (ENS scheme), the discrete formulation of the 
vorticity term provides a global conservation of the enstrophy 
($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent 
flow ($i.e.$ $\chi$=$0$), but does not conserve the total kinetic energy. It is given by:
\begin{equation} \label{Eq_dynvor_ens}
\left\{ 
\begin{aligned}
{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} 
                                & {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i, j+1/2}    \\
{- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}  
                                & {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2, j}  
\end{aligned} 
 \right.
\end{equation} 

%-------------------------------------------------------------
%                 energy conserving scheme
%-------------------------------------------------------------
\subsubsection{Energy conserving scheme (\np{ln\_dynvor\_ene}=true)}
\label{DYN_vor_ene}

The kinetic energy conserving scheme (ENE scheme) conserves the global 
kinetic energy but not the global enstrophy. It is given by:
\begin{equation} \label{Eq_dynvor_ene}
\left\{   \begin{aligned}
{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
                            \;  \overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
{- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
                            \;  \overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} }
\end{aligned}    \right.
\end{equation} 

%-------------------------------------------------------------
%                 mix energy/enstrophy conserving scheme
%-------------------------------------------------------------
\subsubsection{Mixed energy/enstrophy conserving scheme (\np{ln\_dynvor\_mix}=true) }
\label{DYN_vor_mix}

The mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the 
two previous schemes is used. It consists of the ENS scheme (\ref{Eq_dynvor_ens}) 
to the relative vorticity term, and of the ENE scheme (\ref{Eq_dynvor_ene}) applied 
to the planetary vorticity term.
\begin{equation} \label{Eq_dynvor_mix}
\left\{ {     \begin{aligned}
 {+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^{\,i} 
 \; {\overline{\overline {\left( {e_{1v}\,e_{3v}\;v} \right)}} }^{\,i,j+1/2} -\frac{1}{e_{1u} }
 \; {\overline {\left( {\frac{f}{e_{3f} }} \right) 
 \;\overline {\left( {e_{1v}\,e_{3v}\;v} \right)} ^{\,i+1/2}} }^{\,j} } \\
{-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta }{e_{3f} }} \right)} }^j
 \; {\overline{\overline {\left( {e_{2u}\,e_{3u}\;u} \right)}} }^{\,i+1/2,j} +\frac{1}{e_{2v} }
 \; {\overline {\left( {\frac{f}{e_{3f} }} \right)
 \;\overline {\left( {e_{2u}\,e_{3u}\;u} \right)} ^{\,j+1/2}} }^{\,i} } \hfill
\end{aligned}     } \right.
\end{equation} 

%-------------------------------------------------------------
%                 energy and enstrophy conserving scheme
%-------------------------------------------------------------
\subsubsection{Energy and enstrophy conserving scheme (\np{ln\_dynvor\_een}=true) }
\label{DYN_vor_een}

In both the ENS and ENE schemes, it is apparent that the combination of $i$ and $j$ 
averages of the velocity allows for the presence of grid point oscillation structures 
that will be invisible to the operator. These structures are \textit{computational modes} 
that will beat least partly damped by the momentum diffusive operator ($i.e.$ the 
subgrid-scale advection), but not by the resolved advection term. These two schemes
therefore do not contribute to dump grid point noise in the horizontal velocity field, 
which results in more noise in vertical velocity field, an undesired feature. This is a well-known 
characteristics of $C$-grid discretization where $u$ and $v$ are located at different grid point,
a price to pay to avoid a double averaging on the pressure gradient term as in $B$-grid. 
To circumvent this, Adcroft (ADD REF HERE) 
we have proposed to introduce a second velocity ... blahblah.... 

Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....

A very nice solution to that problem was proposed by \citet{Arakawa_Hsu_MWR90}. The idea is
to get rid of the double averaging by considering triad combinations of vorticity. 
It is noteworthy that this solution is conceptually quite similar to the one proposed by
\citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusive operator.

The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified 
for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme. 
Let first provides the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 
\begin{equation} \label{Eq_pot_vor}
q  = \frac{\zeta +f} {e_{3f} }
\end{equation}
where the relative vorticity is defined by (\ref{Eq_divcur_cur}), the Coriolis parameter 
is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 
\begin{equation} \label{Eq_een_e3f}
e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2}
\end{equation}

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!ht] \label{Fig_DYN_een_triad}
\begin{center}
\includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_DYN_een_triad.pdf}
\caption{Triads used in the energy and enstrophy conserving scheme (een) for 
$u$-component (upper panel) and $v$-component (lower panel).}
\end{center}
\end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Note that a key point in \eqref{Eq_een_e3f} is that the averaging in \textbf{i}- and 
\textbf{j}- directions uses the masked vertical scale factor but is always divided by 
$4$, not by the sum of the mask at $T$-point. This preserves the continuity of 
$e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and 
extends by continuity the value of $e_{3f}$ in the land areas. This feature is essential for 
$z$-coordinate with partial step.


Then, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at $T$-point as 
the following triad combinations of the neighbouring potential vorticities defined at f-point 
(Fig.~\ref{Fig_DYN_een_triad}): 
\begin{equation} \label{Q_triads}
_i^j \mathbb{Q}^{i_p}_{j_p}
= \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)
\end{equation}
where the indices $i_p$ and $k_p$ take the following value: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. 

The vorticity terms are represented as: 
\begin{equation} \label{Eq_dynvor_een}
\left\{ {
\begin{aligned}
 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }   \sum_{\substack{i_p,\,k_p}} 
                         {^{i+1/2-i_p}_j}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{1v}\,e_{3v} \;v  \right)^{i+1/2-i_p}_{j+j_p}   \\
 - q\,e_3 \, u     &\equiv -\frac{1}{e_{2v} }    \sum_{\substack{i_p,\,k_p}} 
                         {^i_{j+1/2-j_p}}  \mathbb{Q}^{i_p}_{j_p}  \left( e_{2u}\,e_{3u} \;u  \right)^{i+i_p}_{j+1/2-j_p}   \\
\end{aligned} 
} \right.
\end{equation} 

This EEN scheme in fact combines the conservation properties of ENS and ENE schemes. 
It conserves both total energy and potential enstrophy in the limit of horizontally 
nondivergent flow ($i.e.$ $\chi$=$0$) (see  Appendix~\ref{Apdx_C_vor_zad}). 
Applied to a realistic ocean configuration, it has been shown that its larger mantis
leads to a significant reduction of the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. 
Furthermore, used in combination with partial step representation of bottom topography,
it improves the interaction between current and topography, leading to a larger
topostrophy of the flow  \citep{Barnier_al_OD06, Penduff_al_OS07}. 

%--------------------------------------------------------------------------------------------------------------
%           Kinetic Energy Gradient term
%--------------------------------------------------------------------------------------------------------------
\subsection   [Kinetic Energy Gradient term (\textit{dynkeg})]
         {Kinetic Energy Gradient term (\mdl{dynkeg})}
\label{DYN_keg}

As demonstarted in Appendix~\ref{Apdx_C}, there is a single discrete formulation 
of the kinetic energy gradient term that, together with the formulation chosen for 
the vertical advection (see below), conserves the total kinetic energy:
\begin{equation} \label{Eq_dynkeg}
\left\{ \begin{aligned}
 -\frac{1}{2 \; e_{1u} }  & \ \delta _{i+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]   \\
 -\frac{1}{2 \; e_{2v} }  & \ \delta _{j+1/2} \left[ {\overline {u^2}^{\,i} + \overline{v^2}^{\,j}} \right]    
\end{aligned} \right.
\end{equation} 

%--------------------------------------------------------------------------------------------------------------
%           Vertical advection term
%--------------------------------------------------------------------------------------------------------------
\subsection   [Vertical advection term (\textit{dynzad}) ]
         {Vertical advection term (\mdl{dynzad}) }
\label{DYN_zad}

The discrete formulation of the vertical advection, together with the formulation 
chosen for the gradient of kinetic energy (KE) term, conserves the total kinetic 
energy. Indeed, the change of KE due to the vertical advection is exactly 
balanced by the change of KE due to the gradient of KE (see Appendix~\ref{Apdx_C}).
\begin{equation} \label{Eq_dynzad}
\left\{     \begin{aligned}
-\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}  \\
-\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} 
\end{aligned}         \right.
\end{equation} 

% ================================================================
% Coriolis and Advection : flux form
% ================================================================
\section{Coriolis and Advection: flux form}
\label{DYN_adv_cor_flux}
%------------------------------------------nam_dynadv----------------------------------------------------
\namdisplay{namdyn_adv} 
%-------------------------------------------------------------------------------------------------------------

In the flux form (as in the vector invariant form), the Coriolis and momentum 
advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity 
appearing in their expressions is centred in time (\textit{now} velocity). At the 
lateral boundaries either free slip, no slip or partial slip boundary conditions 
are applied following Chap.\ref{LBC}.


%--------------------------------------------------------------------------------------------------------------
%           Coriolis plus curvature metric terms
%--------------------------------------------------------------------------------------------------------------
\subsection   [Coriolis plus curvature metric terms (\textit{dynvor}) ]
         {Coriolis plus curvature metric terms (\mdl{dynvor}) }
\label{DYN_cor_flux}

In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis 
parameter has been modified to account for the "metric" term. This altered 
Coriolis parameter is thus discretised at $f$-points. It is given by: 
\begin{multline} \label{Eq_dyncor_metric}
f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\
   \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta _{i+1/2} \left[ {e_{2u} } \right]  
                                                                 -  \overline u ^{j+1/2}\delta _{j+1/2} \left[ {e_{1u} } \right]  }  \ \right)
\end{multline} 

Any of the (\ref{Eq_dynvor_ens}), (\ref{Eq_dynvor_ene}) and (\ref{Eq_dynvor_een}) 
schemes can be used to compute the product of the Coriolis parameter and the 
vorticity. However, the energy-conserving scheme  (\ref{Eq_dynvor_een}) has 
exclusively been used to date. This term is evaluated using a leapfrog scheme, 
$i.e.$ the velocity is centred in time (\textit{now} velocity).

%--------------------------------------------------------------------------------------------------------------
%           Flux form Advection term
%--------------------------------------------------------------------------------------------------------------
\subsection   [Flux form Advection term (\textit{dynadv}) ]
         {Flux form Advection term (\mdl{dynadv}) }
\label{DYN_adv_flux}

The discrete expression of the advection term is given by :
\begin{equation} \label{Eq_dynadv}
\left\{ 
\begin{aligned}
\frac{1}{e_{1u}\,e_{2u}\,e_{3u}} 
\left(      \delta _{i+1/2} \left[ \overline{e_{2u}\,e_{3u}\;u }^{i       }  \ u_t      \right]    
          + \delta _{j       } \left[ \overline{e_{1u}\,e_{3u}\;v }^{i+1/2}  \ u_f      \right] \right.  \ \;   \\
\left.   + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{i+1/2}  \ u_{uw} \right] \right)   \\
\\
\frac{1}{e_{1v}\,e_{2v}\,e_{3v}} 
\left(     \delta _{i       } \left[ \overline{e_{2u}\,e_{3u }\;u }^{j+1/2} \ v_f       \right] 
         + \delta _{j+1/2} \left[ \overline{e_{1u}\,e_{3u }\;v }^{i       } \ v_t       \right] \right.  \ \, \, \\
\left.  + \delta _{k      } \left[ \overline{e_{1w}\,e_{2w}\;w}^{j+1/2} \ v_{vw}  \right] \right) \\
\end{aligned}
\right.
\end{equation}

Two advection schemes are available: a $2^{nd}$ order centered finite 
difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS. 
The latter is described in \citet{Shchepetkin_McWilliams_OM05}. The schemes are 
selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}. 
In flux form, the schemes differ by the choice of a space and time interpolation to 
define the value of $u$ and $v$ at the centre of each face of $u$- and $v$-cells, 
$i.e.$ at the $T$-, $f$-, and $uw$-points for $u$ and at the $f$-, $T$- and 
$vw$-points for $v$. 

%-------------------------------------------------------------
%                 2nd order centred scheme
%-------------------------------------------------------------
\subsubsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_dynadv\_cen2}=true)}
\label{DYN_adv_cen2}

In the centered $2^{nd}$ order formulation, the velocity is evaluated as the 
mean of the two neighbouring points :
\begin{equation} \label{Eq_dynadv_cen2}
\left\{     \begin{aligned}
 u_T^{cen2} &=\overline u^{i }       \quad &  u_F^{cen2} &=\overline u^{j+1/2}  \quad &  u_{uw}^{cen2} &=\overline u^{k+1/2}   \\
 v_F^{cen2} &=\overline v ^{i+1/2} \quad & v_F^{cen2} &=\overline v^j      \quad &  v_{vw}^{cen2} &=\overline v ^{k+1/2}  \\
\end{aligned}      \right.
\end{equation} 

The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive 
($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be 
used in conjunction with an explicit diffusion operator to produce a sensible solution. 
The associated time-stepping is performed using a leapfrog scheme in conjunction 
with an Asselin time-filter, so $u$ and $v$ are the \emph{now} velocities.

%-------------------------------------------------------------
%                 UBS scheme
%-------------------------------------------------------------
\subsubsection{Upstream Biased Scheme (UBS) (\np{ln\_dynadv\_ubs}=true)}
\label{DYN_adv_ubs}

The UBS advection scheme is an upstream biased third order scheme based on 
an upstream-biased parabolic interpolation. For example, the evaluation of 
$u_T^{ubs} $ is done as follows:
\begin{equation} \label{Eq_dynadv_ubs}
u_T^{ubs} =\overline u ^i-\;\frac{1}{6}   \begin{cases}
      u"_{i-1/2}&    \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  \geqslant 0$ }    \\
      u"_{i+1/2}&    \text{if $\ \overline{e_{2u}\,e_{3u} \ u}^i  < 0$ }
\end{cases}
\end{equation}
where $u"_{i+1/2} =\delta _{i+1/2} \left[ {\delta _i \left[ u \right]} \right]$. This results 
in a dissipatively dominant ($i.e.$ hyper-diffusive) truncation error \citep{Shchepetkin_McWilliams_OM05}. 
The overall performance of the advection scheme is similar to that reported in 
\citet{Farrow1995}. It is a relatively good compromise between accuracy and 
smoothness. It is not a \emph{positive} scheme, meaning that false extrema are 
permitted. But the amplitudes of the false extrema are significantly reduced over 
those in the centred second order method. As the scheme already includes 
a diffusive component, it can be used without explicit  lateral diffusion on moment 
($i.e.$ \np{ln\_dynldf\_lap}=\np{ln\_dynldf\_bilap}=false), and it is recommended to do so.

The UBS scheme is not used in all directions. In the vertical, the centred $2^{nd}$ 
order evaluation of the advection is preferred, $i.e.$ $u_{uw}^{ubs}$ and 
$u_{vw}^{ubs}$ in \eqref{Eq_dynadv_cen2} are used. UBS is diffusive and is 
associated with vertical mixing of momentum. \gmcomment{ gm  pursue the 
sentence:Since vertical mixing of momentum is a source term of the TKE equation...  }

For stability reasons,  the first term in (\ref{Eq_dynadv_ubs}), which corresponds 
to a second order centred scheme, is evaluated using the \textit{now} velocity 
(centred in time), while the second term, which is the diffusive part of the scheme, 
is evaluated using the \textit{before} velocity (forward in time). This is discussed 
by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.

Note that the UBS and Quadratic Upstream Interpolation for Convective Kinematics 
(QUICK) schemes only differ by one coefficient. Substituting $1/6$ with $1/8$ in 
(\ref{Eq_dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
This option is not available through a namelist parameter, since the $1/6$ coefficient 
is hard coded. Nevertheless it is quite easy to make the substitution in 
\mdl{dynadv\_ubs} module and obtain a QUICK scheme.

Note also that in the current version of \mdl{dynadv\_ubs}, there is also the 
possibility of using a $4^{th}$ order evaluation of the advective velocity as in 
ROMS. This is an error and should be suppressed soon.
%%%
\gmcomment{action :  this have to be done}
%%%

% ================================================================
%           Hydrostatic pressure gradient term
% ================================================================
\section  [Hydrostatic pressure gradient (\textit{dynhpg})]
      {Hydrostatic pressure gradient (\mdl{dynhpg})}
\label{DYN_hpg}
%------------------------------------------nam_dynhpg---------------------------------------------------
\namdisplay{namdyn_hpg} 
%-------------------------------------------------------------------------------------------------------------

The key distinction between the different algorithms used for the hydrostatic 
pressure gradient is the vertical coordinate used, since HPG is a \emph{horizontal} 
pressure gradient, $i.e.$ computed along geopotential surfaces. As a result, any 
tilt of the surface of the computational levels will require a specific treatment to 
compute the hydrostatic pressure gradient.

The hydrostatic pressure gradient term is evaluated either using a leapfrog scheme, 
$i.e.$ the density appearing in its expression is centred in time (\emph{now} $\rho$), or 
a semi-implcit scheme. At the lateral boundaries either free slip, no slip or partial slip 
boundary conditions are applied.

%--------------------------------------------------------------------------------------------------------------
%           z-coordinate with full step
%--------------------------------------------------------------------------------------------------------------
\subsection   [$z$-coordinate with full step (\np{ln\_dynhpg\_zco}) ]
         {$z$-coordinate with full step (\np{ln\_dynhpg\_zco}=true)}
\label{DYN_hpg_zco}

The hydrostatic pressure can be obtained by integrating the hydrostatic equation 
vertically from the surface. However, the pressure is large at great depth while its 
horizontal gradient is several orders of magnitude smaller. This may lead to large 
truncation errors in the pressure gradient terms. Thus, the two horizontal components 
of the hydrostatic pressure gradient are computed directly as follows:

for $k=km$ (surface layer, $jk=1$ in the code)
\begin{equation} \label{Eq_dynhpg_zco_surf}
\left\{ \begin{aligned}
               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k=km} 
&= \frac{1}{2} g \   \left. \delta _{i+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
                  \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k=km} 
&= \frac{1}{2} g \   \left. \delta _{j+1/2} \left[  e_{3w} \ \rho \right] \right|_{k=km}   \\
\end{aligned} \right.
\end{equation} 

for $1<k<km$ (interior layer)
\begin{equation} \label{Eq_dynhpg_zco}
\left\{ \begin{aligned}
               \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k} 
&=             \left. \delta _{i+1/2} \left[  p^h         \right] \right|_{k-1} 
+    \frac{1}{2}\;g\;   \left. \delta _{i+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
                  \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k} 
&=                \left. \delta _{j+1/2} \left[  p^h            \right] \right|_{k-1} 
+    \frac{1}{2}\;g\;   \left. \delta _{j+1/2} \left[  e_{3w} \ \overline {\rho}^{k+1/2} \right] \right|_{k}   \\
\end{aligned} \right.
\end{equation} 

Note that the $1/2$ factor in (\ref{Eq_dynhpg_zco_surf}) is adequate because of 
the definition of $e_{3w}$ as the vertical derivative of the scale factor at the surface 
level ($z=0$). Note also that in case of variable volume level (\key{vvl} defined), the 
surface pressure gradient is included in \eqref{Eq_dynhpg_zco_surf} and \eqref{Eq_dynhpg_zco} 
through the space and time variations of the vertical scale factor $e_{3w}$.

%--------------------------------------------------------------------------------------------------------------
%           z-coordinate with partial step
%--------------------------------------------------------------------------------------------------------------
\subsection   [$z$-coordinate with partial step (\np{ln\_dynhpg\_zps})]
         {$z$-coordinate with partial step (\np{ln\_dynhpg\_zps}=true)}
\label{DYN_hpg_zps}

With partial bottom cells, tracers in horizontally adjacent cells generally live at 
different depths. Before taking horizontal gradients between these tracer points, 
a linear interpolation is used to approximate the deeper tracer as if it actually lived 
at the depth of the shallower tracer point. 

Apart from this modification, the horizontal hydrostatic pressure gradient evaluated 
in the $z$-coordinate with partial step is exactly as in the pure $z$-coordinate case. 
As explained in detail in section \S\ref{TRA_zpshde}, the nonlinearity of pressure 
effects in the equation of state is such that it is better to interpolate temperature and 
salinity vertically before computing the density. Horizontal gradients of temperature 
and salinity are needed for the TRA modules, which is the reason why the horizontal 
gradients of density at the deepest model level are computed in module \mdl{zpsdhe} 
located in the TRA directory and described in \S\ref{TRA_zpshde}.

%--------------------------------------------------------------------------------------------------------------
%           s- and s-z-coordinates
%--------------------------------------------------------------------------------------------------------------
\subsection{$s$- and $z$-$s$-coordinates}
\label{DYN_hpg_sco}

Pressure gradient formulations in $s$-coordinate have been the subject of a vast 
literature ($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}). 
A number of different pressure gradient options are coded, but they are not yet fully 
documented or tested. 

$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true, 
\np{ln\_dynhpg\_hel}=true)
\begin{equation} \label{Eq_dynhpg_sco}
\left\{ \begin{aligned}
 - \frac{1}                   {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  p^h  \right] 
+ \frac{g\; \overline {\rho}^{i+1/2}}  {\rho_o \, e_{1u}} \;   \delta _{i+1/2} \left[  z_t   \right]    \\
 - \frac{1}                   {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  p^h  \right]  
+ \frac{g\; \overline {\rho}^{j+1/2}}  {\rho_o \, e_{2v}} \;   \delta _{j+1/2} \left[  z_t   \right]    \\
\end{aligned} \right.
\end{equation} 

Where the first term is the pressure gradient along coordinates, computed as in 
\eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of 
the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point 
($e_{3w}$). The version \np{ln\_dynhpg\_hel}=true has been added by Aike 
Beckmann and involves a redefinition of the relative position of $T$-points relative 
to $w$-points. 

$\bullet$ Weighted density Jacobian (WDJ) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=true)

$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05} 
(\np{ln\_dynhpg\_djc}=true)

$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true)

Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume 
formulation is activated (\key{vvl}) because in that case, even with a flat bottom, 
the coordinate surfaces are not horizontal but follow the free surface 
\citep{Levier2007}. The other pressure gradient options are not yet available.

%--------------------------------------------------------------------------------------------------------------
%           Time-scheme
%--------------------------------------------------------------------------------------------------------------
\subsection   [Time-scheme (\np{ln\_dynhpg\_imp}) ]
         {Time-scheme (\np{ln\_dynhpg\_imp}= true/false)}
\label{DYN_hpg_imp}

The default time differencing scheme used for the horizontal pressure gradient is 
a leapfrog scheme and therefore the density used in all discrete expressions given 
above is the  \textit{now} density, computed from the \textit{now} temperature and 
salinity. In some specific cases (usually high resolution simulations over an ocean 
domain which includes weakly stratified regions) the physical phenomenum that 
controls the time-step is internal gravity waves (IGWs). A semi-implicit scheme for 
doubling the stability limit associated with IGWs can be used \citep{Brown_Campana_MWR78, 
Maltrud1998}. It involves the evaluation of the hydrostatic pressure gradient as an 
average over the three time levels $t-\rdt$, $t$, and $t+\rdt$ ($i.e.$  
\textit{before},  \textit{now} and  \textit{after} time-steps), rather than at central 
time level $t$ only, as in the standard leapfrog scheme. 

$\bullet$ leapfrog scheme (\np{ln\_dynhpg\_imp}=true):

\begin{equation} \label{Eq_dynhpg_lf}
\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\frac{1}{\rho _o \,e_{1u} }\delta _{i+1/2} \left[ {p_h^t } \right]
\end{equation}

$\bullet$ semi-implicit scheme (\np{ln\_dynhpg\_imp}=true):
\begin{equation} \label{Eq_dynhpg_imp}
\frac{u^{t+\rdt}-u^{t-\rdt}}{2\rdt} = \;\cdots \;
   -\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]
\end{equation}

The semi-implicit time scheme \eqref{Eq_dynhpg_imp} is made possible without 
significant additional computation since the density can be updated to time level 
$t+\rdt$ before computing the horizontal hydrostatic pressure gradient. It can 
be easily shown that the stability limit associated with the hydrostatic pressure 
gradient doubles using \eqref{Eq_dynhpg_imp} compared to that using the 
standard leapfrog scheme \eqref{Eq_dynhpg_lf}. Note that \eqref{Eq_dynhpg_imp} 
is equivalent to applying a time filter to the pressure gradient to eliminate high 
frequency IGWs. Obviously, when using \eqref{Eq_dynhpg_imp}, the doubling of 
the time-step is achievable only if no other factors control the time-step, such as 
the stability limits associated with advection or diffusion.

In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=true. 
In this case, we choose to apply the time filter to temperature and salinity used in 
the equation of state, instead of applying it to the hydrostatic pressure or to the 
density, so that no additional storage array has to be defined. The density used to 
compute the hydrostatic pressure gradient (whatever the formulation) is evaluated 
as follows:
\begin{equation} \label{Eq_rho_flt}
   \rho^t = \rho( \widetilde{T},\widetilde {S},z_t)
 \quad     \text{with}  \quad 
   \widetilde{X} = 1 / 4 \left(  X^{t+\rdt} +2 \,X^t + X^{t-\rdt}  \right)
\end{equation}

Note that in the semi-implicit case, it is necessary to save the filtered density, an 
extra three-dimensional field, in the restart file to restart the model with exact 
reproducibility. This option is controlled by  \np{nn\_dynhpg\_rst}, a namelist parameter.

% ================================================================
% Surface Pressure Gradient
% ================================================================
\section  [Surface pressure gradient (\textit{dynspg}) ]
      {Surface pressure gradient (\mdl{dynspg})}
\label{DYN_spg}
%-----------------------------------------nam_dynspg----------------------------------------------------
\namdisplay{namdyn_spg} 
%------------------------------------------------------------------------------------------------------------

$\ $\newline      %force an empty line

%%%
The 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}.

%%%


The form of the surface pressure gradient term depends on how the user wants to handle 
the fast external gravity waves that are solution of the analytical equation (\S\ref{PE_hor_pg}). 
Three formulation are available, all controlled by a CPP key (ln\_dynspg\_xxx):
a explicit formulation which required a small time step ;
a filtered free surface formulation which allows a larger time step by adding a filtering 
term in the momentum equation ; 
and a plit-explicit free surface formulation, described below, which also allows a larger time step.

The extra term introduced in the filtered method is calculated 
implicitly, so that a solver is used to compute it and that the update of the $next$ 
velocities is done in module \mdl{dynspg\_flt} and not in \mdl{dynnxt}.



%--------------------------------------------------------------------------------------------------------------
% Explicit free surface formulation
%--------------------------------------------------------------------------------------------------------------
\subsection{Explicit free surface (\key{dynspg\_exp})}
\label{DYN_spg_exp}

In the explicit free surface formulation (\key{dynspg\_exp} defined), the model time step 
is chosen to be small enough to describe the external gravity waves (typically a few tens of seconds). 
The surface pressure gradient, evaluated using a leap-frog scheme ($i.e.$ centered in time),
is thus simply given by :
\begin{equation} \label{Eq_dynspg_exp}
\left\{ \begin{aligned}
 - \frac{1}{e_{1u}\,\rho_o} \;   \delta _{i+1/2} \left[  \,\rho \,\eta\,  \right]   \\
 - \frac{1}{e_{2v}\,\rho_o} \;   \delta _{j+1/2} \left[  \,\rho \,\eta\,  \right]  
\end{aligned} \right.
\end{equation} 

Note that in the non-linear free surface case ($i.e.$ \key{vvl} defined), the surface pressure 
gradient is already included in the momentum tendency  through the level thickness variation 
when computing the hydrostatic pressure gradient. Thus, nothing is done in the \mdl{dynspg\_exp} module.

%--------------------------------------------------------------------------------------------------------------
% Split-explict free surface formulation
%--------------------------------------------------------------------------------------------------------------
\subsection{Split-Explicit free surface (\key{dynspg\_ts})}
\label{DYN_spg_ts}

In the split-explicit free surface formulation, also called time-splitting formulation 
(\key{dynspg\_ts} defined)


The split-explicit free surface formulation used in \NEMO follows the one 
proposed by \citet{Griffies_Bk04}. The general idea is to solve the free surface 
equation and the associated barotropic velocity equations with a smaller time 
step than $\rdt$, the time step used for the three dimensional prognostic 
variables (Fig.\ref {Fig_DYN_dynspg_ts}). 
The size of the small time step, $\Delta_e$ (the external mode or barotropic time step)
 is provided through \np{nn\_baro} namelist parameter as: 
$\Delta_e = \Delta / nn\_baro$.
 

%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
\begin{figure}[!t] \label{Fig_DYN_dynspg_ts}
\begin{center}
\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf}
\caption{Schematic of the split-explicit time stepping scheme for the external 
and internal modes. Time increases to the right. 
Internal mode time steps (which are also the model time steps) are denoted 
by $t-\rdt$, $t, t+\rdt$, and $t+2\rdt$. 
The curved line represents a leap-frog time step, and the smaller time 
steps $N \rdt_e=\frac{3}{2}\rdt$ are denoted by the zig-zag line. 
The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ 
represents the interaction between the external and internal motions. 
While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog 
integration carries the external mode variables (surface height and vertically 
integrated velocity) from $t$ to $t+\frac{3}{2} \rdt$ using N external time 
steps of length $\rdt_e$. Time averaging the external fields over the 
$\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated 
velocity and the sea surface height at the model timestep $t+\rdt$. 
These averaged values are used to update \textbf{M}(t) with both the surface 
pressure gradient and the Coriolis force, therefore providing the $t+\rdt$
velocity.  The model time stepping scheme can then be achieved by a baroclinic 
leap-frog time step that carries the surface height from $t-\rdt$ to $t+\rdt$.  }
\end{center}
\end{figure}
%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >

The split-explicit formulation has a damping effect on external gravity waves, 
which is weaker damping than for the filtered free surface but still significant as 
shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave. 

%>>>>>===============
\gmcomment{               %%% copy from griffies Book 

\textbf{title: Time stepping the barotropic system }

Assume knowledge of the full velocity and tracer fields at baroclinic time $\tau$. Hence, 
we can update the surface height and vertically integrated velocity with a leap-frog 
scheme using the small barotropic time step $\rdt$. We have 

\begin{equation} \label{DYN_spg_ts_eta}
\eta^{(b)}(\tau,t_{n+1}) - \eta^{(b)}(\tau,t_{n+1}) (\tau,t_{n-1})
   = 2 \rdt \left[-\nabla \cdot \textbf{U}^{(b)}(\tau,t_n) + \text{EMP}_w(\tau) \right] 
\end{equation}
\begin{multline} \label{DYN_spg_ts_u}
\textbf{U}^{(b)}(\tau,t_{n+1}) - \textbf{U}^{(b)}(\tau,t_{n-1})  \\
   = 2\rdt \left[ - f \textbf{k} \times \textbf{U}^{(b)}(\tau,t_{n}) 
   - H(\tau) \nabla p_s^{(b)}(\tau,t_{n}) +\textbf{M}(\tau) \right]
\end{multline}
\

In 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)}$ 
and $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 
that sets the barotropic time steps via 
\begin{equation} \label{DYN_spg_ts_t}
t_n=\tau+n\rdt   
\end{equation}
with $n$ an integer. The density scaled surface pressure is evaluated via 
\begin{equation} \label{DYN_spg_ts_ps}
p_s^{(b)}(\tau,t_{n}) = \begin{cases}
   g \;\eta_s^{(b)}(\tau,t_{n}) \;\rho(\tau)_{k=1}) / \rho_o  &      \text{non-linear case} \\
   g \;\eta_s^{(b)}(\tau,t_{n})  &      \text{linear case} 
   \end{cases}
\end{equation}
To get started, we assume the following initial conditions 
\begin{equation} \label{DYN_spg_ts_eta}
\begin{split}
\eta^{(b)}(\tau,t_{n=0}) &= \overline{\eta^{(b)}(\tau)}
\\
\eta^{(b)}(\tau,t_{n=1}) &= \eta^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0} 
\end{split}
\end{equation}
with 
\begin{equation} \label{DYN_spg_ts_etaF}
 \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
\end{equation}
the time averaged surface height taken from the previous barotropic cycle. Likewise, 
\begin{equation} \label{DYN_spg_ts_u}
\textbf{U}^{(b)}(\tau,t_{n=0}) = \overline{\textbf{U}^{(b)}(\tau)}   \\
\\
\textbf{U}(\tau,t_{n=1}) = \textbf{U}^{(b)}(\tau,t_{n=0}) + \rdt \ \text{RHS}_{n=0}    
\end{equation}
with 
\begin{equation} \label{DYN_spg_ts_u}
 \overline{\textbf{U}^{(b)}(\tau)} 
   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau-\rdt,t_{n})
\end{equation}
the time averaged vertically integrated transport. Notably, there is no Robert-Asselin time filter used in the barotropic portion of the integration. 

Upon reaching $t_{n=N} = \tau + 2\Delta \tau$ , the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at baroclinic time $\tau + \Delta \tau$ 
\begin{equation} \label{DYN_spg_ts_u}
\textbf{U}(\tau+\rdt) = \overline{\textbf{U}^{(b)}(\tau+\rdt)} 
   = \frac{1}{N+1} \sum\limits_{n=0}^N\textbf{U}^{(b)}(\tau,t_{n})
\end{equation}
The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using the following form 

\begin{equation} \label{DYN_spg_ts_ssh}
\eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]  
\end{equation}

 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). 
 
In general, some form of time filter is needed to maintain integrity of the surface 
height field due to the leap-frog splitting mode in equation \ref{DYN_spg_ts_ssh}. We 
have tried various forms of such filtering, with the following method discussed in 
\cite{Griffies_al_MWR01} chosen due to its stability and reasonably good maintenance of 
tracer conservation properties (see Section ??) 

\begin{equation} \label{DYN_spg_ts_sshf}
\eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)} 
\end{equation}
Another approach tried was 

\begin{equation} \label{DYN_spg_ts_sshf2}
\eta^{F}(\tau-\Delta) = \eta(\tau) 
   + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
                + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
\end{equation}

which is useful since it isolates all the time filtering aspects into the term multiplied 
by $\alpha$. This isolation allows for an easy check that tracer conservation is exact when 
eliminating 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. 

}            %%end gm comment (copy of griffies book)

%>>>>>===============


%--------------------------------------------------------------------------------------------------------------
% Filtered free surface formulation
%--------------------------------------------------------------------------------------------------------------
\subsection{Filtered free surface (\key{dynspg\_flt})}
\label{DYN_spg_fltp}

The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation. 
The extra term introduced in the equations (see {\S}I.2.2) is solved implicitly. 
The elliptic solvers available in the code are documented in \S\ref{MISC}.

%% gm %%======>>>>   given here the discrete eqs provided to the solver

Note that in the linear free surface formulation (\key{vvl} not defined), the ocean depth 
is time-dependent and so is the matrix to be inverted. It is computed once for all the
ocean time steps. 

% ================================================================
% Lateral diffusion term
% ================================================================
\section  [Lateral diffusion term (\textit{dynldf})]
      {Lateral diffusion term (\mdl{dynldf})}
\label{DYN_ldf}
%------------------------------------------nam_dynldf----------------------------------------------------
\namdisplay{namdyn_ldf} 
%-------------------------------------------------------------------------------------------------------------

The options available for lateral diffusion are for the choice of  laplacian 
(rotated or not) or biharmonic operators. The coefficients may be constant 
or spatially variable; the description of the coefficients is found in the chapter 
on lateral physics (Chap.\ref{LDF}). The lateral diffusion of momentum is 
evaluated using a forward scheme, $i.e.$ the velocity appearing in its expression 
is the \textit{before} velocity in time, except for the pure vertical component 
that appears when a tensor of rotation is used. This latter term is solved 
implicitly together with the vertical diffusion term (see \S\ref{DOM_nxt}) 

At the lateral boundaries either free slip, no slip or partial slip boundary 
conditions are applied according to the user's choice (see Chap.\ref{LBC}).

% ================================================================
\subsection   [Iso-level laplacian operator (\np{ln\_dynldf\_lap}) ]
         {Iso-level laplacian operator (\np{ln\_dynldf\_lap}=true)}
\label{DYN_ldf_lap}

For lateral iso-level diffusion, the discrete operator is: 
\begin{equation} \label{Eq_dynldf_lap}
\left\{ \begin{aligned}
 D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta _{i+1/2} \left[ {A_T^{lm} 
\;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta _j \left[ 
{A_f^{lm} \;e_{3f} \zeta } \right] \\ 
\\
 D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta _{j+1/2} \left[ {A_T^{lm} 
\;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta _i \left[ 
{A_f^{lm} \;e_{3f} \zeta } \right] \\ 
\end{aligned} \right.
\end{equation} 

As explained in \S\ref{PE_ldf}, this formulation (as the gradient of a divergence 
and curl of the vorticity) preserves symmetry and ensures a complete 
separation between the vorticity and divergence parts of the momentum diffusion. 
Note that in the full step $z$-coordinate (\key{zco} is defined), $e_{3u} =e_{3v} =e_{3f}$ 
so that they cancel in the rotational part of \eqref{Eq_dynldf_lap}.

%--------------------------------------------------------------------------------------------------------------
%           Rotated laplacian operator
%--------------------------------------------------------------------------------------------------------------
\subsection   [Rotated laplacian operator (\np{ln\_dynldf\_iso}) ]
         {Rotated laplacian operator (\np{ln\_dynldf\_iso}=true)}
\label{DYN_ldf_iso}

A rotation of the lateral momentum diffusive operator is needed in several cases: 
for iso-neutral diffusion in $z$-coordinate (\np{ln\_dynldf\_iso}=true) and for 
either iso-neutral (\np{ln\_dynldf\_iso}=true) or geopotential 
(\np{ln\_dynldf\_hor}=true) diffusion in $s$-coordinate. In the partial step 
case, coordinates are horizontal excepted at the deepest level and no 
rotation is performed when \np{ln\_dynldf\_hor}=true. The diffusive operator 
is defined simply as the divergence of down gradient momentum fluxes on each 
momentum component. It must be emphasized that this formulation ignores 
constraints on the stress tensor such as symmetry. The resulting discrete 
representation is:
\begin{equation} \label{Eq_dyn_ldf_iso}
\begin{split}
 D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} } \\
&  \left\{\quad  {\delta _{i+1/2} \left[ {A_T^{lm}  \left( 
    {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta _{i}[u]
   -e_{2t} \; r_{1t} \,\overline{\overline {\delta _{k+1/2}[u]}}^{\,i,\,k}}
 \right)} \right]}   \right.
\\ 
& \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f} 
}\,\delta _{j+1/2} [u] - e_{1f}\, r_{2f} 
\,\overline{\overline {\delta _{k+1/2} [u]}} ^{\,j+1/2,\,k}} 
\right)} \right] 
\\ 
&\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline 
{\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} } 
\right.} \right. 
\\ 
&  \ \qquad \qquad \qquad \quad\ 
- e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2}
\\ 
& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 
+\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2} 
\right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} 
\\
\\
 D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} }    \\
&  \left\{\quad  {\delta _{i+1/2} \left[ {A_f^{lm}  \left( 
    {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta _{i+1/2}[v]
   -e_{2f} \; r_{1f} \,\overline{\overline {\delta _{k+1/2}[v]}}^{\,i+1/2,\,k}}
 \right)} \right]}   \right.
\\ 
& \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t} 
}\,\delta _{j} [v] - e_{1t}\, r_{2t} 
\,\overline{\overline {\delta _{k+1/2} [v]}} ^{\,j,\,k}} 
\right)} \right] 
\\ 
& \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline 
{\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. 
\\
&  \ \qquad \qquad \qquad \quad\ 
- e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2}
\\ 
& \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\ 
+\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2} 
\right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\} 
 \end{split}
\end{equation}
where $r_1$ and $r_2$ are the slopes between the surface along which the 
diffusive operator acts and the surface of computation ($z$- or $s$-surfaces). 
The way these slopes are evaluated is given in the lateral physics chapter 
(Chap.\ref{LDF}).

%--------------------------------------------------------------------------------------------------------------
%           Iso-level bilaplacian operator
%--------------------------------------------------------------------------------------------------------------
\subsection   [Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap})]
         {Iso-level bilaplacian operator (\np{ln\_dynldf\_bilap}=true)}
\label{DYN_ldf_bilap}

The lateral fourth order operator formulation on momentum is obtained by 
applying \eqref{Eq_dynldf_lap} twice. It requires an additional assumption on 
boundary conditions: the first derivative term normal to the coast depends on 
the free or no-slip lateral boundary conditions chosen, while the third 
derivative terms normal to the coast are set to zero (see Chap.\ref{LBC}).
%%%
\gmcomment{add a remark on the the change in the position of the coefficient}
%%%

% ================================================================
%           Vertical diffusion term
% ================================================================
\section  [Vertical diffusion term (\mdl{dynzdf})]
      {Vertical diffusion term (\mdl{dynzdf})}
\label{DYN_zdf}
%----------------------------------------------namzdf------------------------------------------------------
\namdisplay{namzdf} 
%-------------------------------------------------------------------------------------------------------------

The large vertical diffusion coefficient found in the surface mixed layer together 
with high vertical resolution implies that in the case of explicit time stepping there 
would be too restrictive a constraint on the time step. Two time stepping schemes 
can be used for the vertical diffusion term : $(a)$ a forward time differencing 
scheme (\np{ln\_zdfexp}=true) using a time splitting technique 
(\np{nn\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme 
(\np{ln\_zdfexp}=false) (see \S\ref{DOM_nxt}). Note that namelist variables 
\np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 

The formulation of the vertical subgrid scale physics is the same whatever 
the vertical coordinate is. The vertical diffusion operators given by 
\eqref{Eq_PE_zdf} take the following semi-discrete space form:
\begin{equation} \label{Eq_dynzdf}
\left\{   \begin{aligned}
D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta _k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
                              \ \delta _{k+1/2} [\,u\,]         \right]     \\
\\
D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta _k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
                              \ \delta _{k+1/2} [\,v\,]         \right]
\end{aligned}   \right.
\end{equation} 
where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and 
diffusivity coefficients. The way these coefficients are evaluated 
depends on the vertical physics used (see \S\ref{ZDF}).

The surface boundary condition on momentum is given by the stress exerted by 
the wind. At the surface, the momentum fluxes are prescribed as the boundary 
condition on the vertical turbulent momentum fluxes,
\begin{equation} \label{Eq_dynzdf_sbc}
\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
    = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v }
\end{equation}
where $\left( \tau _u ,\tau _v \right)$ are the two components of the wind stress 
vector in the (\textbf{i},\textbf{j}) coordinate system. The high mixing coefficients 
in the surface mixed layer ensure that the surface wind stress is distributed in 
the vertical over the mixed layer depth. If the vertical mixing coefficient 
is small (when no mixed layer scheme is used) the surface stress enters only 
the top model level, as a body force. The surface wind stress is calculated 
in the surface module routines (SBC, see Chap.\ref{SBC})

The turbulent flux of momentum at the bottom of the ocean is specified through 
a bottom friction parameterisation (see \S\ref{ZDF_bfr})

% ================================================================
% External Forcing
% ================================================================
\section{External Forcings}
\label{DYN_forcing}

Besides the surface and bottom stresses (see the above section) which are 
introduced as boundary conditions on the vertical mixing, two other forcings 
enter the dynamical equations. 

One is the effect of atmospheric pressure on the ocean dynamics (to be 
introduced later). 

Another forcing term is the tidal potential, which will be introduced in the 
reference version soon. 

% ================================================================
% Time evolution term 
% ================================================================
\section  [Time evolution term (\textit{dynnxt})]
      {Time evolution term (\mdl{dynnxt})}
\label{DYN_nxt}

%----------------------------------------------namdom----------------------------------------------------
\namdisplay{namdom} 
%-------------------------------------------------------------------------------------------------------------

The general framework for dynamics time stepping is a leap-frog scheme, 
$i.e.$ a three level centred time scheme associated with an Asselin time filter 
(cf. Chap.\ref{STP}). The scheme is applied to the velocity except when using 
the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in variable 
volume level case (\key{vvl} defined), where it has to be applied to the thickness 
weighted velocity (see \S\ref{Apdx_A_momentum})  

$\bullet$ vector invariant form or linear free surface (\np{ln\_dynhpg\_vec}=true ; not \key{vvl}):
\begin{equation} \label{Eq_dynnxt_vec}
\left\{   \begin{aligned}
&u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t     \\
&u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
\end{aligned}   \right.
\end{equation} 

$\bullet$ flux form and nonlinear free surface (\np{ln\_dynhpg\_vec}=false ; \key{vvl} defined):
\begin{equation} \label{Eq_dynnxt_flux}
\left\{   \begin{aligned}
&\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t     \\
&\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
  +\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]
\end{aligned}   \right.
\end{equation} 
where RHS is the right hand side of the momentum equation, the subscript $f$ 
denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is 
initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = $10^{-3}$.
In both cases, the modified Asselin filter is not applied since a perfect conservation 
is not an issue for momentum equation.

Note that with the filtered free surface, the update of the \textit{next} velocities 
is done in the \mdl{dynsp\_flt} module, and only the swap of arrays 
and Asselin filtering is done in \mdl{dynnxt}.



% ================================================================