source: trunk/NEMO/DOC/BETA/Chapters/Chap_DYN.tex @ 707

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

% add a figure for  dynvor ens, ene latices


$\ $\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 conservative 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 mean 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 in a vorticity part (VOR), a kinetic energy part (KEG) and the vertical advection (ZAD) in the vector invariant formulation, and into Coriolis and advection (COR+ADV) in the flux formulation. The following terms are the pressure gradient contributions (HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient). Contributions from lateral diffusion and vertical diffusion are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules; the latter 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 of the SBC, LDF and ZDF categories and 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 accounting for 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 on the rhs of the 3D momentum equation (\key{trddyn} defined) and of the 2D barotropic vorticity balance (\key{trdvor} defined), as described in Chap.~\ref{MISC}.

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

The vector invariant form of the momentum equations is the most often used in applications of the ocean model. The flux form option (next section) has been introduced recently in version $2$ of NEMO. 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}.

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

Different discretisations of the vorticity term (selected by the \textit{ln\_dynvor\_xxx} namelist variable to true) are available, that conserve potential enstrophy of horizontally non-divergent flow, horizontal kinetic energy, or potential enstrophy for the relative vorticity term and horizontal kinetic energy for the planetary vorticity term (see  \colorbox{yellow} { appendix C } ). The vorticity terms are given below for the general case ($s$-coordinate or partial step topography), but note that in full step $z$-coordinate (\key{zco} defined), $e_{3u} =e_{3v} =e_{3f}$ so that the vertical scale factors disappear.

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

In this case, 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 of 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}=T)}
\label{DYN_vor_ene}

The kinetic energy conserving 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}=T) }
\label{DYN_vor_mix}

In this case, a mixture of the two previous schemes is used. It consists of the enstrophy conserving scheme (\ref{Eq_dynvor_ens}) applied to the relative vorticity term and of the horizontal kinetic energy conserving 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}=T) }
\label{DYN_vor_een}

In this case, the vorticity term is evaluated using the vorticity advection scheme of \citet{Arakawa1990}. This scheme conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow (i.e. $\chi$=0), While this scheme is more complicated that the vorticity advection scheme of \citet{Sadourny1975} and does not conserve potential enstrophy and total energy in general flow, as does a scheme from \citet{Arakawa1981}, it tolerates arbitrarily thin layers. This feature is essential for simulating either outcropping isopycnals or large amplitude topography. 

The \citet{Arakawa1990} vorticity advection scheme for a single layer is modified for spherical coordinates as described by \citet{Arakawa1981}. 

The potential vorticity, defined at $f$-point, is: 
\begin{equation} \label{Eq_pot_vor}
q_f  = \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]{./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. 

The vorticity terms are represented as: 
\begin{equation} \label{Eq_dynvor_een}
\left\{ {
\begin{aligned}
 +q\,e_3 \, v  &\equiv +\frac{1}{e_{1u} }  \left[ 
{{\begin{array}{*{20}c}
      {\,\ \ a_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} } 
   {\,+\,b_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i-1/2}  } \\
 \\
     {  +\,c_{j-1/2}^{i   }  \left( {e_{1v} e_{3v} \ v} \right)_{j    }^{i+1/2}         } 
   {\,+\,d_{j+1/2}^{i   } \left( {e_{1v} e_{3v} \ v} \right)_{j+1}^{i+1/2} } \\
\end{array} }} \right] \\ 
\\
-q\,e_3 \,u       &\equiv -\frac{1}{e_{2v} }  \left[ 
{{\begin{array}{*{20}c}
   {\,\ \ a_{j-1/2}^{i   }  \left( {e_{2u} e_{3v} \ u} \right)_{j+1}^{i+1/2} } 
   {\,+\,b_{j-1/2}^{i+1}  \left( {e_{2u} e_{3v} \ u} \right)_{j+1/2}^{i+1} } \\
 \\
      {  +\,c_{j+1/2}^{i+1} \left( {e_{2u} e_{3v} \ u} \right)_{j+1/2}^{i+1} } 
   {\,+\,d_{j+1/2}^{i   }  \left( {e_{2u} e_{3v} \ u} \right)_{j+1/2}^{i   } } \\
\end{array} }} \right]
\end{aligned} 
} \right.
\end{equation} 
where $a$, $b$, $c$ and $d$ are triad combinations of the neighbouring potential vorticities (Fig. \ref{Fig_DYN_een_triad}): 
\begin{equation} \label{Eq_een_triads}
\left\{ 
\begin{aligned}
 a_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j+1/2}^{i+1} + q_{j+1 /2}^i + q_{j-1/2}^i  \right)    \\ 
 \\
 b_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j+1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i   \right)     \\ 
\\
 c_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j-1/2}^{i-1} +q_{j+1/2}^i +q_{j-1/2}^i   \right)     \\ 
\\
 d_{\,j+1/2}^i & = \frac{1}  {12} \left( q_{j-1/2}^{i+1} +q_{j+1/2}^i +q_{j-1/2}^i  \right)     \\ 
\end{aligned} 
\right.
\end{equation}

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

There is a single discrete formulation of the kinetic energy gradient term, that conserves the total kinetic energy together with the formulation chosen for the vertical advection (see below). 
\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 (\mdl{dynzad}) }
\label{DYN_zad}

The discrete formulation of the vertical advection term conserves the total kinetic energy together with the formulation chosen for the gradient of kinetic energy (KE). 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 \colorbox{yellow}{Annexe 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{nam_dynadv} 
%-------------------------------------------------------------------------------------------------------------

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 (\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 (\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 ($2^{nd}$ order centered finite difference scheme, CEN2, or a $3^{rd}$ order upstream biased scheme, UBS). The latter is described in \citet{Sacha2005}. The schemes are selected using the namelist variable \textit{ln\_dynadv\_xxx} In flux form, the advection 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 and $f$-, $T$- and $vw$-points for $u$ and $v$, respectively. 

%-------------------------------------------------------------
%                 2nd order centred scheme
%-------------------------------------------------------------
\subsubsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_dynadv\_cen2}=T)}
\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}=T)}
\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{Sacha2005}. 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 their amplitude is significantly reduced over the centred second order method. 

The UBS scheme is not used in all the direction. In the vertical, it has been preferred to keep the centred $2^{nd}$ order evaluation of the advection, i.e. $u_{uw}^{ubs}$ and $u_{vw}^{ubs}$ in (\ref{Eq_dynadv_cen2}) are used. Indeed, UBS is diffusive, it is associated with vertical mixing of momentum. Since vertical mixing of momentum is source term of the TKE equation{\ldots}

For stability reasons, in (\ref{Eq_dynadv_ubs}), the first term 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{Webb1998} in the context of the Quick advection scheme.

NB 1: 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{Webb1998}. 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.

NB 2 : In the current version of \mdl{dynadv\_ubs}, there is also a possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. This is an error and should be suppressed soon.

% ================================================================
%           Hydrostatic pressure gradient term
% ================================================================
\section{Hydrostatic pressure gradient (\mdl{dynhpg})}
\label{DYN_hpg}
%------------------------------------------nam_dynhpg---------------------------------------------------
\namdisplay{nam_dynhpg} 
\namdisplay{namflg} 
%-------------------------------------------------------------------------------------------------------------
\colorbox{red} {Suppress the namflg namelist and incorporate it in the namhpg namelist}

The key distinction between the different algorithms is the vertical coordinate used, since HPG is a 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} velocity), 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{\textit{z}-coordinate with full step (\np{ln\_dynhpg\_zco}=T)}
\label{DYN_hpg_zco}

The hydrostatic pressure can be obtained by integrating the hydrostatic equation vertically from the surface. Nevertheless, the pressure is quite large at great depth while its horizontal gradient is several orders of magnitude smaller. This may lead to large truncation error 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)$. 

%--------------------------------------------------------------------------------------------------------------
%           z-coordinate with partial step
%--------------------------------------------------------------------------------------------------------------
\subsection{\textit{z}-coordinate with partial step (\np{ln\_dynhpg\_zps}=T )}
\label{DYN_hpg_zps}

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

It is necessary to do so at the last ocean level; otherwise the horizontal hydrostatic pressure gradient evaluated in z-coordinate with partial step is exactly as in 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 vertically temperature and salinity before computing the density. Horizontal gradients of 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{\textit{s}- and \textit{s-z} coordinates}
\label{DYN_hpg_sco}

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

Traditional coding (see for example \citet{Madec1996}: (\np{ln\_dynhpg\_sco}, \np{ln\_dynhpg\_hel})
\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 (II.3.2), and $z_T$ is the depth of $T$-point evaluated from the sum of the vertical scale factor at W-point ($e_{3w})$. The version \np{ln\_dynhpg\_hel} has been added by Aike Beckmann and involves a redefinition of the relative position of $T$ points relative to $w$ points. 

Weighted density Jacobian (wdj) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=T)

Density Jacobian with cubic polynomial scheme (djc) \citep{Sacha2003} (\np{ln\_dynhpg\_djc}=T)

Rotated axes scheme (rot) \citep{Thiem2006} (\np{ln\_dynhpg\_rot}=T)

Note that expression (\ref{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 surface 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}=T/F)}
\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 a ocean domain including weakly stratified regions) the physical phenomena that control the time-step are internal gravity waves (IGWs). A semi-implicit scheme for doubling the stability limit associated with IGWs can be used \citep{Brown1978, Maltrud1998}. It consists in evaluating the hydrostatic pressure gradient as an average over the three time levels $t-\Delta t$, $t$, and $t+\Delta t$ (i.e.  \textit{before},  \textit{now} and  \textit{after} time-steps), rather than at central time level $t$ only as in standard leapfrog scheme. 

leapfrog scheme (\np{ln\_dynhpg\_imp}=F):

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

semi-implicit scheme (\np{ln\_dynhpg\_imp}=T):
\begin{equation} \label{Eq_dynhpg_imp}
\frac{u^{t+\Delta t}-u^{t-\Delta t}}{2\Delta t}
=\;\cdots \;-\frac{1}{\rho _o \,e_{1u} } \delta _{i+1/2} \left[ \frac{ p_h^{t+\Delta t} +2p_h^t 
+p_h^{t-\Delta t} } { 4 }  \right]
\end{equation}

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

In practice, the semi-implicit scheme is used when \np{ln\_dynhpg\_imp}=T.
 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{\bullet} = \frac{  \bullet^{t+\Delta t} +2 \,\bullet^t + \bullet^{t-\Delta t}  } {4}
\end{equation}

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

% ================================================================
% Surface Pressure Gradient
% ================================================================
\section{Surface pressure gradient (\mdl{dynspg})}
\label{DYN_hpg_spg}
%-----------------------------------------nam_dynspg----------------------------------------------------
\namdisplay{nam_dynspg} 
%------------------------------------------------------------------------------------------------------------
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 or rigid lid) and the variable volume case (nonlinear free surface, \key{vvl} is active). In the linear free surface case (\S\ref{PE_free_surface}) and rigid lid (\S\ref{PE_rigid_lid}), the vertical scale factors $e_{3}$ are fixed in time, while in the nonlinear case (\S\ref{PE_free_surface}) they are time-dependent. 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}.

%--------------------------------------------------------------------------------------------------------------
% Linear free surface formulation
%--------------------------------------------------------------------------------------------------------------
\subsection{Linear free surface formulation (\textbf{key\_exp} or \textbf{\_ts} or \textbf{\_flt})}
\label{DYN_spg_linear}

In the linear free surface formulation, the sea surface height is assumed to be small compared to the thickness of the ocean levels, so that $(a)$ the time evolution of the sea surface height becomes a linear equation, and $(b)$ the thickness of the ocean levels is assumed to be constant with time. As mentioned in (\S\ref{PE_free_surface}) the linearization affects the conservation of tracers. 

%-------------------------------------------------------------
% Explicit
%-------------------------------------------------------------
\subsubsection{Explicit (\key{dynspg\_exp})}
\label{DYN_spg_exp}

In the explicit free surface formulation, the model time step is chosen small enough to describe the external gravity waves (typically a few ten seconds). The sea surface height is given by :
\begin{equation} \label{Eq_dynspg_ssh}
\frac{\partial \eta }{\partial t}\equiv \frac{\text{EMP}}{\rho _w }+\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)} 
\end{equation}

where EMP is the surface freshwater budget (evaporation minus precipitation, and minus river runoffs (if the later are introduced as a surface freshwater flux, see \S\ref{SBC}) expressed in $Kg.m^{-2}.s^{-1}$, and $\rho _w =1,000\,Kg.m^{-3}$ is the volumic mass of pure water. The sea-surface height is evaluated using a leapfrog scheme in combination with an Asselin time filter, i.e. the velocity appearing in (\ref{Eq_dynspg_ssh}) is centred in time (\textit{now} velocity). 

The surface pressure gradient, also evaluated using a leap-frog scheme, is then simply given by :
\begin{equation} \label{Eq_dynspg_exp}
\left\{ \begin{aligned}
 - \frac{1}                      {e_{1u}} \; \delta _{i+1/2} \left[  \,\eta\,  \right]    \\
 \\
 - \frac{1}                      {e_{2v}} \; \delta _{j+1/2} \left[  \,\eta\,  \right]  
\end{aligned} \right.
\end{equation} 

Consistent with the linearization, a $\left. \rho \right|_{k=1} / \rho _o$ factor is omitted in (\ref{Eq_dynspg_exp}). 

%-------------------------------------------------------------
% Split-explicit time-stepping
%-------------------------------------------------------------
\subsubsection{Split-explicit time-stepping (\key{dynspg\_ts})}
\label{DYN_spg_ts}
%--------------------------------------------namdom----------------------------------------------------
\namdisplay{namdom} 
%--------------------------------------------------------------------------------------------------------------
The split-explicit free surface formulation used in OPA follows the one proposed by \citet{Griffies2004}. The general idea is to solve the free surface equation with a small time step, while the three dimensional prognostic variables are solved with a longer time step that is a multiple of \np{rdtbt} (Figure III.3). 

%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
\begin{figure}[!t] \label{Fig_DYN_dynspg_ts}
\begin{center}
\includegraphics[width=0.90\textwidth]{./Figures/Fig_DYN_dynspg_ts.pdf}
\caption{Schematic of the split-explicit time stepping scheme for the barotropic and baroclinic modes, after \citet{Griffies2004}. Time increases to the right. Baroclinic time steps are denoted by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. The curved line represents a leap-frog time step, and the smaller barotropic time steps $N \Delta t=2\Delta t$ are denoted by the zig-zag line. The vertically integrated forcing \textbf{M}(t) computed at baroclinic time step t represents the interaction between the barotropic and baroclinic motions. While keeping the total depth, tracer, and freshwater forcing fields fixed, a leap-frog integration carries the surface height and vertically integrated velocity from t to $t+2 \Delta t$ using N barotropic time steps of length $\Delta t$. Time averaging the barotropic fields over the N+1 time steps (endpoints included) centers the vertically integrated velocity at the baroclinic timestep $t+\Delta t$. A baroclinic leap-frog time step carries the surface height to $t+\Delta t$ using the convergence of the time averaged vertically integrated velocity taken from baroclinic time step t. }
\end{center}
\end{figure}
%>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >

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

%-------------------------------------------------------------
% Filtered formulation 
%-------------------------------------------------------------
\subsubsection{Filtered formulation (\key{dynspg\_flt})}
\label{DYN_spg_flt}

The filtered formulation follows the \citet{Roullet2000} 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}. The amplitude of the extra term is given by the namelist variable \np{rnu}. The default value is 1, as recommended by \citet{Roullet2000}

\colorbox{red}{\np{rnu}=1 to be suppressed from namelist !}

%-------------------------------------------------------------
% Non-linear free surface formulation 
%-------------------------------------------------------------
\subsection{Non-linear free surface formulation (\key{vvl})}
\label{DYN_spg_vvl}

In the non-linear free surface formulation, the variations of volume are fully taken into account. This option is presented in a report \citep{Levier2007} available on the NEMO web site. The three time-stepping methods (explicit, split-explicit and filtered) are the same as in \S\ref{DYN_spg_linear} except that the ocean depth is now time-dependent. In particular, this means that in filtered case, the matrix to be inverted has to be recomputed at each time-step.

%--------------------------------------------------------------------------------------------------------------
%           Rigid-lid formulation
%--------------------------------------------------------------------------------------------------------------
\subsection{Rigid-lid formulation (\key{dynspg\_rl})}
\label{DYN_spg_rl}

With the rigid lid formulation, an elliptic equation has to be solved for 
the barotropic streamfunction. For consistency this equation is obtained by 
taking the discrete curl of the discrete vertical sum of the discrete 
momentum equation: 
\begin{equation}\label{Eq_dynspg_rl}
\frac{1}{\rho _o }\nabla _h p_s \equiv \left( {{\begin{array}{*{20}c}
 {\overline M_u +\frac{1}{H\;e_2 } \delta_ j \left[ \partial_t \psi \right]}      \\
 \\
 {\overline M_v -\frac{1}{H\;e_1 }  \delta_i  \left[ \partial_t \psi \right]}        \\
\end{array} }} \right)
\end{equation}

Here ${\rm {\bf M}}= \left( M_u,M_v \right)$ represents the collected contributions of nonlinear, 
viscous and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn} and the overbar indicates a vertical average over the whole water column (i.e. from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid). The time derivative of $\psi$ is the solution of an elliptic equation:
\begin{multline} \label{Eq_bsf}
   \delta_{i+1/2} \left[ \frac{e_{2v}}{H_v\;e_{1v}} \delta_{i} \left[  \partial_t \psi \right] \right]
+ \delta_{j+1/2} \left[ \frac{e_{1u}}{H_u\;e_{2u}} \delta_{j} \left[  \partial_t \psi \right] \right]
\\ =
  \delta_{i+1/2} \left[ e_{2v} M_v  \right]
- \delta_{j+1/2} \left[ e_{1u} M_u  \right]
\end{multline}

The elliptic solvers available in the code are documented in 
\S\ref{MISC}). The boundary conditions must be given on all separate 
landmasses (islands), which is done by integrating the vorticity equation 
around each island. This requires identifying each island in the bathymetry 
file, a cumbersome task. This explains why the rigid lid option is not 
recommended with complex domains such as the global ocean. Parameters jpisl 
(number of islands) and jpnisl (maximum number of points per island) of the 
par\_oce.h90 file are related to this option.


% ================================================================
% Lateral diffusion term
% ================================================================
\section{Lateral diffusion term (\mdl{dynldf})}
\label{DYN_ldf}
%------------------------------------------nam_dynldf----------------------------------------------------
\namdisplay{nam_dynldf} 
%-------------------------------------------------------------------------------------------------------------
The options available for lateral diffusion are 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 following user's choice (see Chap. \S\ref{LBC}).

% ================================================================
\subsection{Iso-level laplacian operator (\np{ln\_dynldf\_lap}=T)}
\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. Note that in pure 
z-coordinate (\key{zco} defined), $e_{3u} =e_{3v} =e_{3f}$ so that they disappear 
from the rotational part of \eqref{Eq_dynldf_lap}.

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

A rotation of lateral momentum diffusive operator is needed for isoneutral 
diffusion in $z$-coordinates (\np{ln\_dynldf\_iso}=T) and for either 
isoneutral (\np{ln\_dynldf\_iso}=T) or geopotential 
(\np{ln\_dynldf\_hor}=T) diffusion in $s$-coordinates. In the partial step 
case, coordinates are horizontal excepted at the deepest level and no 
rotation is performed when \np{ln\_dynldf\_hor}=T. The diffusive operator 
is defined simply as the divergence of down gradient momentum fluxes on each 
momentum components. 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 lateral physics chapter (Chap.~\S\ref{LDF}).

%--------------------------------------------------------------------------------------------------------------
%           Iso-level bilaplacian operator
%--------------------------------------------------------------------------------------------------------------
\subsection{Iso-level bilaplacian operator}
\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}).

\colorbox{yellow}{add a remark on the the change in the position of the coefficient}

% ================================================================
%           Vertical diffusion term
% ================================================================
\section{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 a too restrictive constraint on the time step in a pure explicit time stepping 
case. Two time stepping schemes can be used for the vertical diffusion term : $(a)$ a forward time differencing scheme (\np{ln\_zdfexp}=T) using a time splitting technique (\np{n\_zdfexp} $>$ 1) or $(b)$ a backward (or implicit) time differencing scheme (\np{ln\_zdfexp}=F) (see \S\ref{DOM_nxt}). Note that namelist variables \np{ln\_zdfexp} and \np{n\_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 diffusive 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 can be 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.~\S\ref{SBC})

The turbulent flux of momentum at the bottom is specified through a bottom 
friction parameterization (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 (\mdl{dynnxt})}
\label{DYN_nxt}

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

The general framework of dynamics time stepping is a leap-frog scheme, i.e. 
a three level centred time scheme associated with a Asselin time filter (cf. \S\ref{DOM_nxt})
\begin{equation} \label{Eq_dynnxt}
\begin{split}
&u^{t+\Delta t} = u^{t-\Delta t} + 2 \, \Delta t  \ \text{RHS}_u^t   \\
\\
&u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\Delta t} -2u^t+u^{t+\Delta t}} \right]
\end{split}
\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{atfp} (namelist parameter). Its default value is \np{atfp}=0.1.

Note that whith the filtered free surface, the update of the next velocities 
is done in the \mdl{dynsp\_flt} module, and nly the swap of arrays 
and Asselin filter is done in \mdl{dynnxt.}

% ================================================================
% Diagnostic variables
% ================================================================
\section{Diagnostic variables ($\zeta$, $\chi$, $w$)}
\label{DYN_divcur_wzv}

%--------------------------------------------------------------------------------------------------------------
%           z-coordinate with partial step
%--------------------------------------------------------------------------------------------------------------
\subsection{horizontal divergence and relative vorticity (\mdl{divcur})}
\label{DYN_divcur}

The vorticity is defined at 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 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 z-coordinate with full step (\key{zco} defined), $e_{3u} =e_{3v} =e_{3f}$ so that they disappear from \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 respectively for the nonlinear advection and the computation of the vertical velocity. 

%--------------------------------------------------------------------------------------------------------------
%           Vertical Velocity
%--------------------------------------------------------------------------------------------------------------
\subsection{Vertical velocity (\mdl{wzvmod})}
\label{DYN_wzv}

The vertical velocity is computed by an upward integration of the horizontal 
divergence from the bottom :

\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  
\end{aligned}   \right.
\end{equation} 

With a free surface, the top vertical velocity is non-zero, due the 
freshwater forcing and the variations of the free surface elevation. When 
the free surface is linear or with a rigid lid, the upper boundary condition 
applies at a fixed level $z=0$. Note that in the rigid-lid case (\key{dynspg\_rl} defined), the surface 
boundary condition ($\left. w \right|_\text{surface} =0)$ is automatically 
achieved at least at the computer accuracy, due to the discrete expression 
of the surface pressure gradient (\colorbox{yellow}{Appendix C}).

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. 

With the variable volume option, the calculation of the vertical velocity is 
modified (see \citet{Levier2007}, report available on the NEMO web site). 

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