chap_TRA.tex

\[
  % \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 into
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 through the \nam{tra_dmp}{tra\_dmp} namelist variables.
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 \nam{tsd}{tsd} namelist as well as
\np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures are correctly set
(\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 during
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
  artificial walls of a limited domain basin and
\item the computation of the velocity field associated with a given $T$-$S$ field
  (for example 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 that have artificial walls instead of open boundaries.
In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days)
whereas it is zero in the interior of 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 whilst
having 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 within the ocean.
It also has undesirable effects on the ocean convection.
It tends to prevent deep convection and subsequent deep-water formation,
by stabilising the water column too much.

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

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

%% =================================================================================================
\section[Tracer time evolution (\textit{tranxt.F90})]{Tracer time evolution (\protect\mdl{tranxt})}
\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
(\np[=.true.]{ln_linssh}{ln\_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$.

When the vertical mixing is solved implicitly,
the update of the \textit{next} tracer fields is done in \mdl{trazdf} module.
In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} 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[Horizontal derivative in \textit{zps}-coordinate (\textit{zpshde.F90})]{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})}
\label{sec:TRA_zpshde}

\cmtgm{STEVEN: to be consistent with earlier discussion of differencing and averaging operators,
I've changed "derivative" to "difference" and "mean" to "average"}

With partial cells (\np[=.true.]{ln_zps}{ln\_zps}) at bottom and top
(\np[=.true.]{ln_isfcav}{ln\_isfcav}),
in general, tracers in horizontally adjacent cells live at different depths.
Horizontal gradients of tracers are needed for horizontal diffusion
(\mdl{traldf} module) and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
The partial cell properties at the top (\np[=.true.]{ln_isfcav}{ln\_isfcav}) are computed in
the same way as for the bottom.
So, only the bottom interpolation is explained below.

Before taking horizontal gradients between the tracers next to the bottom,
a linear interpolation in the vertical is used to approximate the deeper tracer as if
it actually lived at the depth of the shallower tracer point (\autoref{fig:TRA_Partial_step_scheme}).
For example, for temperature in the $i$-direction the needed interpolated temperature,
$\widetilde T$, is:

\begin{figure}
  \centering
  \includegraphics[width=0.33\textwidth]{TRA_partial_step_scheme}
  \caption[Discretisation of the horizontal difference and average of tracers in
  the $z$-partial step coordinate]{
    Discretisation of the horizontal difference and average of tracers in
    the $z$-partial step coordinate (\protect\np[=.true.]{ln_zps}{ln\_zps}) in
    the case $(e3w_k^{i + 1} - e3w_k^i) > 0$.
    A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$,
    the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points.
    The horizontal difference is then given by:
    $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and the average by:
    $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$.}
  \label{fig:TRA_Partial_step_scheme}
\end{figure}

\[
  \widetilde T = \lt\{
    \begin{alignedat}{2}
      &T^{\, i + 1} &-\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt) }{ e_{3w}^{i + 1} } \; \delta_k T^{i + 1}
      & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\
      &T^{\, i}     &+\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt )}{e_{3w}^i       } \; \delta_k T^{i + 1}
      & \quad \text{if $e_{3w}^{i + 1} <    e_{3w}^i$}
    \end{alignedat}
  \rt.
\]
and the resulting forms for the horizontal difference and the horizontal average value of
$T$ at a $U$-point are:
\begin{equation}
  \label{eq:TRA_zps_hde}
  \begin{split}
    \delta_{i + 1/2} T       &=
    \begin{cases}
      \widetilde T - T^i          & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\
      T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} <    e_{3w}^i
    \end{cases} \\
    \overline T^{\, i + 1/2} &=
    \begin{cases}
      (\widetilde T - T^{\, i}    ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\
      (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} <   e_{3w}^i
    \end{cases}
  \end{split}
\end{equation}

The computation of horizontal derivative of tracers as well as of density is performed once for all at
each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed.
It has to be emphasized that the procedure used to compute the interpolated density,
$\widetilde \rho$, is not the same as that used for $T$ and $S$.
Instead of forming a linear approximation of density,
we compute $\widetilde \rho$ from the interpolated values of $T$ and $S$,
and the pressure at a $u$-point
(in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}):
\[
  % \label{eq:TRA_zps_hde_rho}
  \widetilde \rho = \rho (\widetilde T,\widetilde S,z_u) \quad \text{where~} z_u = \min \lt( z_T^{i + 1},z_T^i \rt)
\]

This is a much better approximation as the variation of $\rho$ with depth (and thus pressure)
is highly non-linear with a true equation of state and thus is badly approximated with
a linear interpolation.
This approximation is used to compute both the horizontal pressure gradient (\autoref{sec:DYN_hpg})
and the slopes of neutral surfaces (\autoref{sec:LDF_slp}).

Note that in almost all the advection schemes presented in this chapter,
both averaging and differencing operators appear.
Yet \autoref{eq:TRA_zps_hde} has not been used in these schemes:
in contrast to diffusion and pressure gradient computations,
no correction for partial steps is applied for advection.
The main motivation is to preserve the domain averaged mean variance of the advected field when
using the $2^{nd}$ order centred scheme.
Sensitivity of the advection schemes to the way horizontal averages are performed in
the vicinity of partial cells should be further investigated in the near future.
\cmtgm{gm :   this last remark has to be done}

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

\end{document}