chap_DYN.tex

\documentclass[../main/NEMO_manual]{subfiles}

\begin{document}

\chapter{Ocean Dynamics (DYN)}
\label{chap:DYN}

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{\textwidth}{l||X|X}
    Release & Author(s) & Modifications \\
    \hline
    {\em   4.0} & {\em ...} & {\em ...} \\
    {\em   3.6} & {\em ...} & {\em ...} \\
    {\em   3.4} & {\em ...} & {\em ...} \\
    {\em <=3.4} & {\em ...} & {\em ...}
  \end{tabularx}
}

\clearpage

Using the representation described in \autoref{chap: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 an average or difference operator is used, the resulting field is multiplied by a mask.

The prognostic ocean dynamics equation can be summarized as follows:
\[
  \text{NXT} = \dbinom	{\text{VOR} + \text{KEG} + \text {ZAD} }
  {\text{COR} + \text{ADV}                       }
  + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
\]
NXT stands for next, referring to the time-stepping.
The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that
are decomposed into either 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
\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.

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

The different options available to the user are managed by namelist variables.
For term \textit{ttt} in the momentum equations, 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 \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 and outputting each tendency term from the 3D momentum equations
(\texttt{trddyn?} defined), as described in \autoref{chap:MISC}.
Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined)
can be derived from the 3D terms.
\cmtgm{STEVEN: not quite sure I've got the sense of the last sentence.
  Does MISC correspond to "extracting tendency terms" or "vorticity balance"?}

%% =================================================================================================
\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)}
\label{sec:DYN_divcur_wzv}

%% =================================================================================================
\subsection[Horizontal divergence and relative vorticity (\textit{divhor.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
\label{subsec:DYN_divcur}

The vorticity is defined at an $f$-point (\ie\ corner point) as follows (used in ):
\begin{equation}
  \label{eq:DYN_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:
\[
  % \label{eq:DYN_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)
\]

Note that although the vorticity has the same discrete expression in $z$- and $s$-coordinates,
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 not necessarily 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.

%% =================================================================================================
\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
\label{subsec:DYN_sshwzv}

The sea surface height is given by:
\begin{equation}
  \label{eq:DYN_spg_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,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation).
If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then
\textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.
In the presence of ice-shelf (\np[=.true.]{ln_isf}{ln\_isf}) 
\textit{emp} can be written as the evaporation minus precipitation, minus the ice-shelf water flux.
In the leapfrog case, the sea-surface height is evaluated using exactly the same time stepping scheme as
the tracer equation \autoref{eq:TRA_nxt}:
a leapfrog scheme in combination with an Asselin time filter,
\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity).
This is of paramount importance.
Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
the sea surface height equation otherwise tracer content will not be conserved
\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}.

The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
taking into account the change of the thickness of the levels:
\begin{equation}
  \label{eq:DYN_wzv}
  \left\{
    \begin{aligned}
      &\left. w \right|_{k_b-1/2} \quad= 0    \qquad \text{where } k_b \text{ is the level just above the sea floor }  	\\
      &\left. w \right|_{k+1/2}     = \left. w \right|_{k-1/2}  +  \left. e_{3t} \right|_{k}\;  \left. \chi \right|_k
      - \frac{1} {2 \rdt} \left(  \left. e_{3t}^{t+1}\right|_{k} - \left. e_{3t}^{t-1}\right|_{k}\right)
    \end{aligned}
  \right.
\end{equation}

In the case of a non-linear free surface (\key{qco}), the top vertical velocity is $-\textit{emp}/\rho_w$,
as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses,
re-orientated downward.
\cmtgm{not sure of this...  to be modified with the change in emp setting}
In the case of a linear free surface(\key{linssh}), the time derivative in \autoref{eq:DYN_wzv} disappears.
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
(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}).

Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates,
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-orientated downwards in the \fortran\ code compared to
the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv}
(see \autoref{subsec:DOM_Num_Index_vertical}).

When \np[=.true.]{ln_zad_Aimp}{ln\_zad\_Aimp},
a proportion of the vertical advection can be treated implicitly (see \autoref{sec:DYN_zdf}) depending on the Courant number. 
This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.

%% =================================================================================================
\section{Coriolis and advection: vector invariant form}
\label{sec:DYN_adv_cor_vect}

\begin{listing}
  \nlst{namdyn_adv}
  \caption{\forcode{&namdyn_adv}}
  \label{lst:namdyn_adv}
\end{listing}

The vector invariant form of the momentum equations is the one most often used in
applications of the \NEMO\ ocean model.
The flux form option (see next section) has been present since version $2$.
Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and
momentum advection terms are evaluated either using a leapfrog scheme or a RK3 scheme.
In the leapfrog case the velocity appearing in these expressions is centred in time (\textit{now} velocity).
In the RK3 case the velocity appearing in these expressions is forward in time (\textit{before} velocity) at stage 1, 
it is is centred in time (\textit{now} velocity) at stage 2 and 3. 
At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following
\autoref{chap:LBC}.

%% =================================================================================================
\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})}
\label{subsec:DYN_vor}

\begin{listing}
  \nlst{namdyn_vor}
  \caption{\forcode{&namdyn_vor}}
  \label{lst:namdyn_vor}