chap_D2D.tex

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

\begin{document}

\chapter{Sea surface height and 2D external mode (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

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

In the present chapter we describe the equations used to compute the the sea surface
height. Recall from \autoref{subsec:MB_free_surface} that allowing a free surface permits
the existence of External Gravity Waves (EGWs). Resolving these fast moving waves explicitly
imposes a severe limit on the model timestep to avoid breaching the CFL condition.
Alternatively, EGWs can be filtered by discretisation of the temporal derivatives using a
split-explicit scheme.  This chapter describes these options controlled via the
\nam{dyn_spg}{dyn\_spg} namelist.

%% =================================================================================================
\section[Sea surface height  (\textit{sshwzv.F90}/\textit{stp2d.F90})]
{Sea surface height (\protect\mdl{sshwzv}/\protect\mdl{stp2d})}
\label{sec:D2D_ssh}

The sea surface height is given by the vertical average of the kinematic surface condition 
(\autoref{subsec:MB_free_surface}):
\begin{equation}
  \label{eq:D2D_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 (i.e. evaporation minus precipitation
+ optional contributions (see \autoref{subsec:D2D_forcing})), expressed in Kg/m$^2$/s (which is
equal to mm/s), and $\rho_w$=1,026~Kg/m$^3$ is the reference density of sea water
(Boussinesq approximation).  The sea-surface height is evaluated using exactly the same
time stepping scheme as the tracer \autoref{eq:TRA_nxt}.  E.g., in the case of
MLF this is a leapfrog scheme in combination with an Asselin time filter, \ie\ the
velocity appearing in \autoref{eq:D2D_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}.

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

With MLF time-stepping there is an explicit free surface formulation
(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}), for which the model time step must be
chosen to be small enough to resolve the external gravity waves (typically a few tens of
seconds).  The sea surface height is evaluated from \autoref{eq:D2D_spg_ssh}
using a leap-frog scheme (\ie\ centered in time) with ssh at the now timestep coming from
the previous time step computation.

%, 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 available with RK3 time-stepping but, if implemented, ssh 
would be computed at each stage from \autoref{eq:D2D_spg_ssh} using values at $n$,
$n+{1\over3}$ and $n+{1\over2}$ respectively.

%% =================================================================================================
\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[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}), 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 
\autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.  

For 'single 1$^{\mathrm st}$', RK3 time-stepping, barotropic velocities and sea surface
height are calculated for the \textit{after} timestep from the \textit{before} values
prior to the first stage of the main RK3 timestepping. This occurs in the new module
\forcode{stp2d.F90} which computes the RHS forcing terms from \textit{before} fields and
then calls the \forcode{dyn_spg_ts} routine to perform the time-stepping as detailed
below. In the MLF time-stepping case, RHS forcing terms are computed inside
\forcode{dyn_spg_ts} itself.

For clarity, in RK3 the RHS forcing terms on momentum are vertically averaged 3D trends of:
\[
  \text{HPG} + \text{LDF} + \text {(COR + RVO)} + \text{KEG} + \text{ZAD} 
\]
for the Vector invariant form (for which the 3D fields are also the RHS terms for
the 1$^{\mathrm st}$ stage RK3 time stepping) and:
\[
  \text{HPG} + \text{LDF} + \text {(COR + MET)} + \text{ADV}
\]
for the flux form. In the latter case ADV must be recomputed at the 1$^{\mathrm st}$ stage.
The new terms in these summaries are: RVO - Relative Vorticity and MET - Metric term.
External forcing is also applied in the form of baroclinic bottom drag and surface winds.
Other external forcings may be optionally applied as detailed in
\autoref{subsec:D2D_forcing}.

For ssh, the external forcing is the net column-average freshwater flux. This is
evaporation minus prepcitation plus optional contributions from other sources as listed in
\autoref{subsec:D2D_forcing}.

%% =================================================================================================
\subsection{Split-explicit time-stepping (\protect\mdl{dynspg\_ts})}
\label{subsec:D2D_spg_ts_ts}

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 with MLF.
    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.

%% =================================================================================================
\subsection{External forcings}
\label{subsec:D2D_forcing}

\noindent The net column-average freshwater flux applied to the ssh equation can also be modified by
the following options:

(1) When \np[=.true.]{ln_rnf}{ln\_rnf} (see \autoref{sec:SBC_rnf}),
river runoff is taken into account when computing the net freshwater flux.

(2) When \np[=.true.]{ln_isf}{ln\_isf} (see \autoref{sec:SBC_isf}), explicit or
parameterised contributions from ice-shelf cavities are taken into account when computing
the net freshwater flux.

(3) When \np[=.true.]{ln_sdw}{ln\_sdw} (see \autoref{subsec:SBC_wave_sdw}), the contribution from
divergence due to Stoke's drift is taken into account when computing the net freshwater
flux.

(4) When \np[=.true.]{ln_isf .AND. ln_isfcpl}{ln\_isf .AND. ln\_isfcpl} (see
\autoref{sec:SBC_isf}), the contribution from coupled ice-sheets is taken into account when
computing the net freshwater flux.

(5) When \np[=.true.]{lk_asminc .AND. ln_sshinc .AND. ln_asmiau}{lk\_asminc .AND.
ln\_sshinc .AND. ln\_asmiau} (see \autoref{sec:ASM_IAU}), the contribution from IAU
weighted ssh increments is taken into account when computing the net freshwater flux.

\vskip 0.5cm

\noindent Besides the surface and bottom baroclinic stresses four other forcings may enter the barotropic
momentum equations.

(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.

(4) When \np[=.true.]{ln_wave .AND. ln_bern_srfc}{ln\_wave .AND. ln\_bern\_srfc} (see
\autoref{subsec:SBC_wave_bhd}) , a depth-uniform wave-induced kinematic pressure term (the
Bernoulli head) is added to the mean pressure.

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

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

\end{document}