chap_D2D.tex

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

\begin{document}

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

\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[Surface pressure gradient (\textit{dynspg.F90})]{Surface pressure gradient (\protect\mdl{dynspg})}
\label{sec:D2D_spg}

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

Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.  The surface
pressure gradient term is related to the representation of the free surface
(\autoref{sec:MB_hor_pg}).  The main distinction is between the fixed volume case (linear
free surface, \ie\ with the inclusion of \key{linssh}) and the variable volume case
(nonlinear free surface, \key{qco}).  In the linear free surface case
(\autoref{subsec:MB_free_surface}) the vertical scale factors $e_{3}$ are fixed in time,
while they are time-dependent in the nonlinear case (\autoref{subsec:MB_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
(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}).  To allow a longer time step for the
three-dimensional equations, one can use a split-explicit free surface
(\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}). In that case, only a quasi-linear
 form of 2d barotropic equations is substepped with a small time increment.

%% =================================================================================================
\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]{Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
\label{subsec:D2D_spg_exp}

In the explicit free surface formulation (\np{ln_dynspg_exp}{ln\_dynspg\_exp} set to true),
the model time step is chosen to be small enough to resolve the external gravity waves
(typically a few tens of seconds).
The surface pressure gradient, evaluated using a leap-frog scheme (\ie\ centered in time),
is thus simply given by :
\begin{equation}
  \label{eq:DYN_spg_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 this option is not yet compatible with RK3 time-stepping.

%% =================================================================================================
\subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})}
\label{subsec:D2D_spg_ts}
%\nlst{namsplit}

The split-explicit free surface formulation used in \NEMO\
(\np{ln_dynspg_ts}{ln\_dynspg\_ts} set to true), also called the time-splitting
formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}.  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 (\autoref{fig:D2D_spg_ts}).  The size of the small time
step, $\rdt_e$ (the external mode or barotropic time step) is provided through the
\np{nn_e}{nn\_e} namelist parameter as: $\rdt_e = \rdt / nn\_e$.  This parameter can be
optionally defined automatically (\np[=.true.]{ln_bt_auto}{ln\_bt\_auto}) considering that
the stability of the barotropic system is essentially controled by external waves
propagation.  Maximum Courant number is in that case time independent, and easily computed
online from the input bathymetry.  Therefore, $\rdt_e$ is adjusted so that the Maximum
allowed Courant number is smaller than \np{rn_bt_cmax}{rn\_bt\_cmax}.

The barotropic mode solves the following equations:
% \begin{subequations}
%  \label{eq:D2D_BT}
\begin{equation}
  \label{eq:D2D_BT_dyn}
  \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}=
  -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h}
  -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}}
\end{equation}
\[
  % \label{eq:D2D_BT_ssh}
  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E
\]
% \end{subequations}
where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing
coupling term between modes, surface atmospheric forcing as well as slowly varying
barotropic terms not explicitly computed to gain efficiency.  The third term on the right
hand side of \autoref{eq:D2D_BT_dyn} represents the bottom stress (see section
\autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.  Temporal
discretization of the system above follows a three-time step Generalized Forward Backward
algorithm detailed in \citet{shchepetkin.mcwilliams_OM05}.  AB3-AM4 coefficients used in
\NEMO\ follow the second-order accurate, "multi-purpose" stability compromise as defined
in \citet{shchepetkin.mcwilliams_ibk09} (see their figure 12, lower left).

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{DYN_dynspg_ts}
  \caption[Split-explicit time stepping scheme for the external and internal modes]{
    Schematic of the split-explicit time stepping scheme for the external and internal modes.
    Time increases to the right.
    In this particular exemple,
    a boxcar averaging window over \np{nn_e}{nn\_e} barotropic time steps is used
    (\np[=1]{nn_bt_flt}{nn\_bt\_flt}) and \np[=5]{nn_e}{nn\_e}.
    Internal mode time steps (which are also the model time steps) are denoted by
    $t-\rdt$, $t$ and $t+\rdt$.
    Variables with $k$ superscript refer to instantaneous barotropic variables,
    $< >$ and $<< >>$ operator refer to time filtered variables using respectively primary
    (red vertical bars) and secondary weights (blue vertical bars).
    The former are used to obtain time filtered quantities at $t+\rdt$ while
    the latter are used to obtain time averaged transports to advect tracers.
    a) Forward time integration:
    \protect\np[=1]{nn_bt_flt}{nn\_bt\_flt}.
    b) Centred time integration:
    \protect\np[=2]{nn_bt_flt}{nn\_bt\_flt}.
    c) Forward time integration with no time filtering (POM-like scheme):
    \protect\np[=3]{nn_bt_flt}{nn\_bt\_flt}.}
  \label{fig:D2D_spg_ts}
\end{figure}

In the default case (\protect\np[=1]{nn_bt_flt}{nn\_bt\_flt}), the external mode is
integrated between \textit{now} and \textit{after} baroclinic time-steps
(\autoref{fig:D2D_spg_ts}a).  To avoid aliasing of fast barotropic motions into three
dimensional equations, time filtering is eventually applied on barotropic quantities.  In
that case, the integration is extended slightly beyond \textit{after} time step to provide
time filtered quantities.  These are used for the subsequent initialization of the
barotropic mode in the following baroclinic step.  Since external mode equations written
at baroclinic time steps finally follow a forward time stepping scheme, asselin filtering
is not applied to barotropic quantities.\\ Alternatively, one can choose to integrate
barotropic equations starting from \textit{before} time step
(\protect\np[=2]{nn_bt_flt}{nn\_bt\_flt}).  Although more computationaly expensive (
\np{nn_e}{nn\_e} additional iterations are indeed necessary), the baroclinic to barotropic
forcing term given at \textit{now} time step become centred in the middle of the
integration window.  It can easily be shown that this property removes part of splitting
errors between modes, which increases the overall numerical robustness.
%references to Patrick Marsaleix' work here. Also work done by SHOM group.


As far as tracer conservation is concerned, barotropic velocities used to advect tracers
must also be updated at \textit{now} time step.  This implies to change the traditional
order of computations in \NEMO: most of momentum trends (including the barotropic mode
calculation) updated first, tracers' after.  Advective barotropic velocities are obtained
by using a secondary set of filtering weights, uniquely defined from the filter
coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}).
Consistency between the time averaged continuity equation and the time stepping of tracers
is here the key to obtain exact conservation.


One can eventually choose to feedback instantaneous values by not using any time filter
(\protect\np[=3]{nn_bt_flt}{nn\_bt\_flt}).
In that case, external mode equations are continuous in time,
\ie\ they are not re-initialized when starting a new sub-stepping sequence.
This is the method used in the POM model for example, the stability being maintained by
refreshing at (almost) each barotropic time step advection and horizontal diffusion terms.
Since the latter terms have not been added in \NEMO\ for computational efficiency,
removing time filtering would be inevitably unstable. One can however add some dissipation, but in the time domain, by slightly modifying the barotropic time stepping coefficients (\citet{demange_JCP19}). This is implemented here through an additional parameter (\np{rn_bt_alpha}{rn\_bt\_alpha}), which controls the amount of temporal diffusion.


\cmtgm{               %%% 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

\[
  % \label{eq: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]
\]
\begin{multline*}
  % \label{eq: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
\[
  % \label{eq:DYN_spg_ts_t}
  t_n=\tau+n\rdt
\]
with $n$ an integer.
The density scaled surface pressure is evaluated via
\[
  % \label{eq: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}
\]
To get started, we assume the following initial conditions
\[
  % \label{eq: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}
\]
with
\[
  % \label{eq:DYN_spg_ts_etaF}
  \overline{\eta^{(b)}(\tau)} = \frac{1}{N+1} \sum\limits_{n=0}^N \eta^{(b)}(\tau-\rdt,t_{n})
\]
the time averaged surface height taken from the previous barotropic cycle.
Likewise,
\[
  % \label{eq: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}
\]
with
\[
  % \label{eq: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})
\]
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\rdt \tau$ ,
the vertically integrated velocity is time averaged to produce the updated vertically integrated velocity at
baroclinic time $\tau + \rdt \tau$
\[
  % \label{eq: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})
\]
The surface height on the new baroclinic time step is then determined via a baroclinic leap-frog using
the following form

\begin{equation}
  \label{eq: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 \autoref{eq:DYN_spg_ts_ssh}.
We have tried various forms of such filtering,
with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
its stability and reasonably good maintenance of tracer conservation properties (see ??).

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

\[
  % \label{eq: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]
\]

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 ?? for more complete discussion).
However, in the general case with a non-zero $\alpha$,
the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.

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

%% =================================================================================================
\section{External forcings}
\label{sec:D2D_forcing}

Besides the surface and bottom stresses (see the above section)
which are introduced as boundary conditions on the vertical mixing,
three other forcings may enter the dynamical equations by affecting the surface pressure gradient.

(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
the atmospheric pressure is taken into account when computing the surface pressure gradient.

(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see \autoref{sec:SBC_TDE}),
the tidal potential is taken into account when computing the surface pressure gradient.

(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and SI3 is used
(\ie\ when the sea-ice is embedded in the ocean),
the snow-ice mass is taken into account when computing the surface pressure gradient.

\cmtgm{ missing : the lateral boundary condition !!!   another external forcing
 }

\subinc{\input{../../global/epilogue}}

\end{document}