chap_TRA.tex

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

\begin{figure}
  \centering
  \includegraphics[width=0.66\textwidth]{TRA_geoth}
  \caption[Geothermal heat flux]{
    Geothermal heat flux (in mW m$^{-2}$) estimated by \cite{lucazeau_GGG19}.}
  \label{fig:TRA_geothermal}
\end{figure}

Usually, it is assumed that there is no heat or salt exchange across the ocean floor,
\ie\ a no-flux boundary condition is applied to active tracers on the ocean floor.
This is the default option in \NEMO, which is implemented using the masking technique.
Nevertheless, there is a non-zero heat flux from the seabed associated with solid earth cooling.
This flux is small compared with surface fluxes (the global mean value is thought to be close to 0.1 W m$^{-2}$ \citep{stein.stein_N92}),
but it continuously warms the ocean and acts on the densest water masses. Including this flux
in a global ocean model can increase the deepest overturning cell
(\ie\ the one associated with the Antarctic Bottom Water) by
a few Sverdrups \citep{emile-geay.madec_OS09}. It should be recalled that including geothermal heating also implies a persistent heat flux out of the surface ocean (and ultimately through the top-of-atmosphere), even at steady state. The net heat flux out of the ocean induced by geothermal forcing is focused in the Southern Ocean, where Antarctic Bottom Water is ultimately upwelled.

The options for this heat exchange at the bottom of the ocean are defined in section \nam{bbc}{bbc}.
The presence of geothermal heating is controlled by
setting the namelist parameter \np{ln_trabbc}{ln\_trabbc} to true.
When \np{nn_geoflx}{nn\_geoflx} is set to 1, a constant geothermal heating is applied,
the value of which is given by the namelist parameter \np{rn_geoflx_cst}{rn\_geoflx\_cst}.
When \np{nn_geoflx}{nn\_geoflx} is set to 2,
a spatially varying geothermal heat flux is applied,
provided by a NetCDF file whose name is defined in the namelist section \nam{bbc}{bbc}. It is recommended to use the state-of-the-art estimate of geothermal heat fluxes by \cite{lucazeau_GGG19}, available on NEMO global grids at https://doi.org/10.17882/103233. This estimate is illustrated in 
\autoref{fig:TRA_geothermal}; the global mean ocean value is 85 mW m$^{-2}$ and the globally integrated flux is 31 TW.

%% =================================================================================================
\section[Bottom boundary layer (\textit{trabbl.F90} - \forcode{ln_trabbl})]{Bottom boundary layer (\protect\mdl{trabbl} - \protect\np{ln_trabbl}{ln\_trabbl})}
\label{sec:TRA_bbl}

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

Options are defined through the \nam{bbl}{bbl} namelist variables.
In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps.
This is not adequate to represent gravity driven downslope flows.
Such flows arise either downstream of sills such as the Strait of Gibraltar or Denmark Strait,
where dense water formed in marginal seas flows into a basin filled with less dense water,
or along the continental slope when dense water masses are formed on a continental shelf.
The amount of entrainment that occurs in these gravity plumes is critical in
determining the density and volume flux of the densest waters of the ocean,
such as Antarctic Bottom Water, or North Atlantic Deep Water.
$z$-coordinate models tend to overestimate the entrainment,
because the gravity flow is mixed vertically by convection as
it goes ''downstairs'' following the step topography,
sometimes over a thickness much larger than the thickness of the observed gravity plume.
A similar problem occurs in the $s$-coordinate when
the thickness of the bottom level varies rapidly downstream of a sill
\citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved.

The idea of the bottom boundary layer (BBL) parameterisation, first introduced by
\citet{beckmann.doscher_JPO97},
is to allow a direct communication between two adjacent bottom cells at different levels,
whenever the densest water is located above the less dense water.
The communication can be by a diffusive flux (diffusive BBL),
an advective flux (advective BBL), or both.
In the current implementation of the BBL, only the tracers are modified, not the velocities.
Furthermore, it only connects ocean bottom cells,
and therefore does not include all the improvements introduced by \citet{campin.goosse_T99}.

%% =================================================================================================
\subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf=1})]{Diffusive bottom boundary layer (\protect\np[=1]{nn_bbl_ldf}{nn\_bbl\_ldf})}
\label{subsec:TRA_bbl_diff}

When applying sigma-diffusion
(\np[=.true.]{ln_trabbl}{ln\_trabbl} and \np{nn_bbl_ldf}{nn\_bbl\_ldf} set to 1),
the diffusive flux between two adjacent cells at the ocean floor is given by
\[
  % \label{eq:TRA_bbl_diff}
  \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T
\]
with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells,
and $A_l^\sigma$ the lateral diffusivity in the BBL.
Following \citet{beckmann.doscher_JPO97}, the latter is prescribed with a spatial dependence,
\ie\ in the conditional form
\begin{equation}
  \label{eq:TRA_bbl_coef}
  A_l^\sigma (i,j,t) =
      \begin{cases}
        A_{bbl} & \text{if~} \nabla_\sigma \rho \cdot \nabla H < 0 \\
        0      & \text{otherwise}
      \end{cases}
\end{equation}
where $A_{bbl}$ is the BBL diffusivity coefficient,
given by the namelist parameter \np{rn_ahtbbl}{rn\_ahtbbl} and
usually set to a value much larger than the one used for lateral mixing in the open ocean.
The constraint in \autoref{eq:TRA_bbl_coef} implies that sigma-like diffusion only occurs when
the density above the sea floor, at the top of the slope, is larger than in the deeper ocean
(see green arrow in \autoref{fig:TRA_bbl}).
In practice, this constraint is applied separately in the two horizontal directions,
and the density gradient in \autoref{eq:TRA_bbl_coef} is evaluated with the log gradient formulation:
\[
  % \label{eq:TRA_bbl_Drho}
  \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S
\]
where $\rho$, $\alpha$ and $\beta$ are functions of
$\overline T^\sigma$, $\overline S^\sigma$ and $\overline H^\sigma$,
the along bottom mean temperature, salinity and depth, respectively.

%% =================================================================================================
\subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv=1,2})]{Advective bottom boundary layer (\protect\np[=1,2]{nn_bbl_adv}{nn\_bbl\_adv})}
\label{subsec:TRA_bbl_adv}

%\sgacomment{
%  "downsloping flow" has been replaced by "downslope flow" in the following
%  if this is not what is meant then "downwards sloping flow" is also a possibility"
%}

\begin{figure}
  \centering
  \includegraphics[width=0.33\textwidth]{TRA_BBL_adv}
  \caption[Advective/diffusive bottom boundary layer]{
    Advective/diffusive Bottom Boundary Layer.
    The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$.
    Red arrows indicate the additional overturning circulation due to the advective BBL.
    The transport of the downslope flow is defined either
    as the transport of the bottom ocean cell (black arrow),
    or as a function of the along slope density gradient.
    The green arrow indicates the diffusive BBL flux directly connecting
    $kup$ and $kdwn$ ocean bottom cells.}
  \label{fig:TRA_bbl}
\end{figure}

%!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
%!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
%!!        i.e. transport proportional to the along-slope density gradient

\cmtgm{This section has to be really written}

When applying an advective BBL (\np[=1..2]{nn_bbl_adv}{nn\_bbl\_adv}),
an overturning circulation is added which connects two adjacent bottom grid-points only if
dense water overlies less dense water on the slope.
The density difference causes dense water to move down the slope.

\begin{description}
\item [{\np[=1]{nn_bbl_adv}{nn\_bbl\_adv}}] the downslope velocity is chosen to
  be the Eulerian ocean velocity just above the topographic step
  (see black arrow in \autoref{fig:TRA_bbl}) \citep{beckmann.doscher_JPO97}.
  It is a \textit{conditional advection}, that is,
  advection is allowed only if dense water overlies less dense water on the slope
  (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$) and if the velocity is directed towards greater depth
  (\ie\ $\vect U \cdot \nabla H > 0$).
\item [{\np[=2]{nn_bbl_adv}{nn\_bbl\_adv}}] the downslope velocity is chosen to be proportional to
  $\Delta \rho$, the density difference between the higher cell and lower cell densities
  \citep{campin.goosse_T99}.
  The advection is allowed only  if dense water overlies less dense water on the slope
  (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$).
  For example, the resulting transport of the downslope flow, here in the $i$-direction
  (\autoref{fig:TRA_bbl}), is simply given by the following expression:
  \[
    % \label{eq:TRA_bbl_Utr}
    u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn})
  \]
  where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as
  \np{rn_gambbl}{rn\_gambbl}, a namelist parameter, and
  \textit{kup} and \textit{kdwn} are the vertical index of the higher and lower cells, respectively.
  The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity,
  and because no direct estimation of this parameter is available, a uniform value has been assumed.
  The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}.
\end{description}

Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using
the upwind scheme.
Such a diffusive advective scheme has been chosen to mimic the entrainment between
the downslope plume and the surrounding water at intermediate depths.
The entrainment is replaced by the vertical mixing implicit in the advection scheme.
Let us consider as an example the case displayed in \autoref{fig:TRA_bbl} where
the density at level $(i,kup)$ is larger than the one at level $(i,kdwn)$.
The advective BBL scheme modifies the tracer time tendency of
the ocean cells near the topographic step by the downslope flow \autoref{eq:TRA_bbl_dw},
the horizontal \autoref{eq:TRA_bbl_hor} and the upward \autoref{eq:TRA_bbl_up} return flows as follows:
\begin{alignat}{5}
  \label{eq:TRA_bbl_dw}
  \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\
  \label{eq:TRA_bbl_hor}
  \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\
  \shortintertext{and for $k =kdw-1,\;..., \; kup$ :}
  \label{eq:TRA_bbl_up}
  \partial_t T^{do}_{k}   &\equiv \partial_t S^{do}_{k}   &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   &&\lt( T^{do}_{k +1} - T^{sh}_{k}   \rt)
\end{alignat}
where $b_t$ is the $T$-cell volume.

Note that the BBL transport, $(u^{tr}_{bbl},v^{tr}_{bbl})$, is available in the model outputs.
It has to be used to compute the effective velocity as well as the effective overturning circulation.

%% =================================================================================================
\section[Tracer damping (\textit{tradmp.F90})]{Tracer damping (\protect\mdl{tradmp})}
\label{sec:TRA_dmp}

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

In some applications it can be useful to add a Newtonian damping term to
the temperature and salinity equations :
\begin{equation}
  \label{eq:TRA_dmp}
    \pd[T]{t} = \cdots - \gamma (T - T_o) \qquad \pd[S]{t} = \cdots - \gamma (S - S_o)
\end{equation}
where $\gamma$ is the inverse of a time scale,
and $T_o$ and $S_o$ are given temperature and salinity fields (usually a climatology).
Options are defined in the namelist section \nam{tra_dmp}{tra\_dmp}.
The restoring term is added when the namelist parameter \np{ln_tradmp}{ln\_tradmp} is set to true.
It also requires that both \np{ln_tsd_init}{ln\_tsd\_init} and
\np{ln_tsd_dmp}{ln\_tsd\_dmp} are set to true in the namelist section\nam{tsd}{tsd} and that
\np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures are correctly defined
(\ie\ that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread},
see \autoref{subsec:SBC_fldread}).
The restoring coefficient $\gamma$ is a three-dimensional array read in
the \rou{tra\_dmp\_init} routine.
The file name is specified by the namelist variable \np{cn_resto}{cn\_resto}.
The \texttt{DMP\_TOOLS} are provided to allow users to generate the netcdf file.

The two main cases in which \autoref{eq:TRA_dmp} is used are
\begin{enumerate*}[label=(\textit{\alph*})]
\item the specification of the boundary conditions along
  the artificial walls of a limited-domain basin, and
\item the computation of the velocity field associated with a given $T$-$S$ field
  (e.g. to build the initial state of a prognostic simulation,
  or to use the resulting velocity field for a passive tracer study).
\end{enumerate*}
The first case applies to regional models with artificial walls instead of open boundaries.
In the vicinity of these walls, $\gamma$ takes on large values (equivalent to a time scale of a few days),
whereas it is zero inside the model domain.
The second case corresponds to the use of the robust diagnostic method \citep{sarmiento.bryan_JGR82}.
It allows us to find the velocity field consistent with the model dynamics, while
keeping a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$).

The robust diagnostic method is very efficient in preventing temperature drift in
intermediate waters, but it produces artificial sources of heat and salt in the ocean.
It also has undesirable effects on the ocean convection.
It tends to prevent deep convection and the subsequent formation of deep-waters
by over-stabilising the water column.

The namelist parameter \np{nn_zdmp}{nn\_zdmp} determines whether damping should be applied
throughout the water column or only below the mixed layer (defined either by a density or $S_o$ criterion).
It is common to set the damping to zero in the mixed layer, as the adjustment time scale there is short
\citep{madec.delecluse.ea_JPO96}.

For generating \textit{resto.nc},
see the documentation for the DMP tools supplied with the source code under \path{./tools/DMP_TOOLS}.

%% =================================================================================================
\section[Tracer time evolution (\textit{traatf.F90})]{Tracer time evolution (\protect\mdl{traatf})}
\label{sec:TRA_nxt}

Options are defined through the \nam{dom}{dom} namelist variables.
The general framework for tracer time stepping is a modified leap-frog scheme
\citep{leclair.madec_OM09}, \ie\ a three level centred time scheme associated with
a Asselin time filter (cf. \autoref{sec:TD_mLF}):
\begin{equation}
  \label{eq:TRA_nxt}
  \begin{alignedat}{5}
    &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\
    &(e_{3t}T)_f^t        &&= (e_{3t}T)^t            &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\
    &                     &&                         &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt]
  \end{alignedat}
\end{equation}
where RHS is the right hand side of the temperature equation,
the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient,
and $S$ is the total forcing applied on $T$ (\ie\ fluxes plus content in mass exchanges).
$\gamma$ is initialized as \np{rn_atfp}{rn\_atfp}, its default value is \forcode{10.e-3}.
Note that the forcing correction term in the filter is not applied in linear free surface
(\key{linssh}) (see \autoref{subsec:TRA_sbc}).
Not also that in constant volume case, the time stepping is performed on $T$,
not on its content, $e_{3t}T$.

Vertical mixing is always solved implicitly, 
with updates to the \textit{next} tracer fields handled in the \mdl{trazdf} module. 
The array swapping and Asselin filtering are performed in the \mdl{traatf} module.

In order to prepare for the computation of the \textit{next} time step,
a swap of tracer arrays is performed: $T^{t - \rdt} = T^t$ and $T^t = T_f$.

%% =================================================================================================
\section[Equation of state (\textit{eosbn2.F90})]{Equation of state (\protect\mdl{eosbn2})}
\label{sec:TRA_eosbn2}

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

%% =================================================================================================
\subsection[Equation of seawater (\forcode{ln_}\{\forcode{teos10,eos80,seos}\})]{Equation of seawater (\protect\np{ln_teos10}{ln\_teos10}, \protect\np{ln_teos80}{ln\_teos80}, or \protect\np{ln_seos}{ln\_seos})}
\label{subsec:TRA_eos}

The \textbf{E}quation \textbf{O}f \textbf{S}eawater (EOS) is
an empirical nonlinear thermodynamic relationship linking
seawater density, $\rho$, to a number of state variables,
most typically temperature, salinity and pressure.
Because density gradients control the pressure gradient force through the hydrostatic balance,
the equation of state provides a fundamental bridge between
the distribution of active tracers and the fluid dynamics.
Nonlinearities of the EOS are of major importance, in particular influencing the circulation through
determination of the static stability below the mixed layer,
thus controlling rates of exchange between the atmosphere and the ocean interior
\citep{roquet.madec.ea_JPO15}.
Therefore an accurate EOS based on either the 1980 equation of state
(EOS-80, \cite{fofonoff.millard_bk83}) or TEOS-10 \citep{ioc.iapso_bk10} standards should
be used anytime a simulation of the real ocean circulation is attempted \citep{roquet.madec.ea_JPO15}.
The use of TEOS-10 is highly recommended because
\begin{enumerate*}[label=(\textit{\roman*})]
\item it is the new official EOS,
\item it is more accurate, being based on an updated database of laboratory measurements, and
\item it uses Conservative Temperature and Absolute Salinity
  (instead of potential temperature and practical salinity for EOS-80),
  both variables being more suitable for use as model variables
  \citep{ioc.iapso_bk10, graham.mcdougall_JPO13}.
\end{enumerate*}
EOS-80 is an obsolescent feature of the \NEMO\ system, kept only for backward compatibility.
For process studies, it is often convenient to use an approximation of the EOS.
To that purposed, a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is also available.

In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed,
with $\rho_o$ a reference density.
Called \textit{rho0} in the code,
$\rho_o$ is set in \mdl{phycst} to a value of \texttt{1,026} $Kg/m^3$.
This is a sensible choice for the reference density used in a Boussinesq ocean climate model,
as, with the exception of only a small percentage of the ocean,
density in the World Ocean varies by no more than 2\% from that value \citep{gill_bk82}.

Options which control the EOS used are defined through the \nam{eos}{eos} namelist variables.

\begin{description}
\item [{\np[=.true.]{ln_teos10}{ln\_teos10}}] the polyTEOS10-bsq equation of seawater
  \citep{roquet.madec.ea_OM15} is used.
  The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
  but it is optimized for a Boussinesq fluid and
  the polynomial expressions have simpler and more computationally efficient expressions for
  their derived quantities which make them more adapted for use in ocean models.
  Note that a slightly higher precision polynomial form is now used
  replacement of the TEOS-10 rational function approximation for hydrographic data analysis
  \citep{ioc.iapso_bk10}.
  A key point is that conservative state variables are used:
  Absolute Salinity (unit: $g/kg$, notation: $S_A$) and
  Conservative Temperature (unit: $\deg{C}$, notation: $\Theta$).
  The pressure in decibars is approximated by the depth in meters.
  With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant.
  It is set to $C_p$ = 3991.86795711963 $J.Kg^{-1}.\deg{K}^{-1}$,
  according to \citet{ioc.iapso_bk10}.
  Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$.
  In particular, the initial state defined by the user have to be given as
  \textit{Conservative} Temperature and \textit{Absolute} Salinity.
  In addition, when using TEOS10, the Conservative SST is converted to potential SST prior to
  either computing the air-sea and ice-sea fluxes (forced mode) or
  sending the SST field to the atmosphere (coupled mode).
\item [{\np[=.true.]{ln_eos80}{ln\_eos80}}] the polyEOS80-bsq equation of seawater is used.
  It takes the same polynomial form as the polyTEOS10,
  but the coefficients have been optimized to accurately fit EOS80 (Roquet, personal comm.).
  The state variables used in both the EOS80 and the ocean model are:
  the Practical Salinity (unit: $psu$, notation: $S_p$) and
  Potential Temperature (unit: $\deg{C}$, notation: $\theta$).
  The pressure in decibars is approximated by the depth in meters.
  With EOS, the specific heat capacity of sea water, $C_p$, is a function of
  temperature, salinity and pressure \citep{fofonoff.millard_bk83}.
  Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which
  is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value.
\item [{\np[=.true.]{ln_seos}{ln\_seos}}] a simplified EOS (S-EOS) inspired by
  \citet{vallis_bk06} is chosen,
  the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.)
  (see also \citet{roquet.madec.ea_JPO15}).
  It provides a simplistic linear representation of both cabbeling and thermobaricity effects which
  is enough for a proper treatment of the EOS in theoretical studies \citep{roquet.madec.ea_JPO15}.
  With such an equation of state there is no longer a distinction between \textit{conservative} and
  \textit{potential} temperature, as well as between \textit{absolute} and
  \textit{practical} salinity.
  S-EOS takes the following expression:
  \begin{gather*}
    % \label{eq:TRA_S-EOS}
    d_a(T,S,z) = \frac{1}{\rho_o} \big[ - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \big.
                                        + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a
                                  \big. - \nu \;                           T_a                  S_a \big] \\
    \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3
  \end{gather*}
  where the computer name of the coefficients as well as their standard value are given in
  \autoref{tab:TRA_SEOS}.
  In fact, when choosing S-EOS, various approximation of EOS can be specified simply by
  changing the associated coefficients.
  Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$
  remove thermobaric effect from S-EOS.
  Setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$
  remove   cabbeling effect from S-EOS.
  Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
\end{description}

\begin{table}
  \centering
  \begin{tabular}{|l|l|l|l|}
    \hline
    coeff.      & computer name                & S-EOS            & description                     \\
    \hline
    $a_0      $ & \np{rn_a0}{rn\_a0}           & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\
    \hline
    $b_0      $ & \np{rn_b0}{rn\_b0}           & $7.6554~10^{-1}$ & linear haline  expansion coeff. \\
    \hline
    $\lambda_1$ & \np{rn_lambda1}{rn\_lambda1} & $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$       \\
    \hline
    $\lambda_2$ & \np{rn_lambda2}{rn\_lambda2} & $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$       \\
    \hline
    $\nu      $ & \np{rn_nu}{rn\_nu}           & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$    \\
    \hline
    $\mu_1    $ & \np{rn_mu1}{rn\_mu1}         & $1.4970~10^{-4}$ & thermobaric coeff. in T         \\
    \hline
    $\mu_2    $ & \np{rn_mu2}{rn\_mu2}         & $1.1090~10^{-5}$ & thermobaric coeff. in S         \\
    \hline
  \end{tabular}
  \caption{Standard value of S-EOS coefficients}
  \label{tab:TRA_SEOS}
\end{table}

%% =================================================================================================
\subsection[Brunt-V\"{a}is\"{a}l\"{a} frequency]{Brunt-V\"{a}is\"{a}l\"{a} frequency}
\label{subsec:TRA_bn2}

An accurate computation of the ocean stability (i.e. of $N$, the Brunt-V\"{a}is\"{a}l\"{a} frequency) is of paramount importance as determine the ocean stratification and
is used in several ocean parameterisations
(namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion,
non-penetrative convection, tidal mixing  parameterisation, iso-neutral diffusion).
In particular, $N^2$ has to be computed at the local pressure
(pressure in decibar being approximated by the depth in meters).
The expression for $N^2$  is given by:
\[
  % \label{eq:TRA_bn2}
  N^2 = \frac{g}{e_{3w}} \lt( \beta \; \delta_{k + 1/2}[S] - \alpha \; \delta_{k + 1/2}[T] \rt)
\]
where $(T,S) = (\Theta,S_A)$ for TEOS10, $(\theta,S_p)$ for TEOS-80, or $(T,S)$ for S-EOS, and,
$\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
The coefficients are a polynomial function of temperature, salinity and depth which
expression depends on the chosen EOS.
They are computed through \textit{eos\_rab}, a \fortran\ function that can be found in \mdl{eosbn2}.

%% =================================================================================================
\subsection{Freezing point of seawater}
\label{subsec:TRA_fzp}

The freezing point of seawater is a function of salinity and pressure \citep{fofonoff.millard_bk83}:
\begin{equation}
  \label{eq:TRA_eos_fzp}
  \begin{gathered}
    T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\
    \text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \text{and~} d = -7.53~10^{-3}
    \end{gathered}
\end{equation}

\autoref{eq:TRA_eos_fzp} is only used to compute the potential freezing point of sea water
(\ie\ referenced to the surface $p = 0$),
thus the pressure dependent terms in \autoref{eq:TRA_eos_fzp} (last term) have been dropped.
The freezing point is computed through \textit{eos\_fzp},
a \fortran\ function that can be found in \mdl{eosbn2}.

%% =================================================================================================
%\subsection{Potential Energy anomalies}
%\label{subsec:TRA_bn2}

%    =====>>>>> TO BE written

%% =================================================================================================
\section[Wave induced transport]{Wave induced transport}
\label{sec:TRA_stoke_drift}

The Stokes drift is a wave-driven mechanism that results in the net transport of mass and momentum, 
defined as the difference between the motion of a fluid parcel (Lagrangian velocity) and 
the flow observed at a fixed point (Eulerian velocity).
Due to the asymmetry in the orbital paths of water particles, a net forward motion occurs. 
Further details are available in \autoref{subsec:SBC_wave_sdw}.

Incorporating Stokes drift is critical for improving ocean circulation models. 
Notably, the tracer advection equation is modified to allow Eulerian ocean models to account for 
unresolved wave effects. 

When simulating waves \np{ln_wave}{ln\_wave} and activating Stokes drift effect \np{ln_sdw}{ln\_sdw} 
three-dimensional Stokes velocity is merely added to the tracers advective transports by \rou{tra\_adv\_trp}, 
a \fortran\ function that can be found in \mdl{traadv}.
Since horizontal velocities are modified, the vertical velocity requires to be recomputed. 

The divergence of the wave tracer flux equals the mean tracer advection induced by the three-dimensional Stokes velocity, 
ensuring the consistence between the continuity equation and tracers evolution equations. 
Thus the barotropic divergence requires to take the stoke drift divergence into account.

Note that Stoke drift velocity also contributes to the shear computation,
the Coriolis term (see \autoref{subsec:SBC_wave_stcor}).
It is also used to compute the vortex force added to the relative vorticity term in the 
Vector Invariant Formulation of the momentum equations (see \autoref{subsec:SBC_wave_vf}).

%% =================================================================================================
\section[Internal wave filtering]{Internal wave filtering}
\label{sec:TRA_shuman}

Internal gravity waves can sometimes disrupt an experiment by making the model unstable, 
restricting the usable time step, or enhancing undesirable vertical mixing. 
They can also introduce complex frequencies that complicate result analysis. 

The studies by \citep{brown.campana_MWR78} and Shuman (1971) demonstrate that 
applying a time-averaging technique to the pressure term can mitigate the impact of these waves, 
thereby allowing the maximum time step to be reached without causing instabilities linked to 
internal gravity waves.
It is equivalent to time-average the pressure term in the momentum right hand side equations and 
the time integrated transports.
In \NEMO\ , for sake of simplicity, we retain the second option.
When \np[=.true.]{ln_shuman}{ln\_shuman} transports are time-averaged in \rou{tra\_adv\_trp}, 
a \fortran\ function that can be found in \mdl{traadv}.
Since horizontal velocities are modified, the vertical velocity requires to be recomputed.

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

\end{document}