chap_ZDF.tex

%% =================================================================================================
\subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf_}\{\forcode{tke,gls,osm}\})]{Handling convection with turbulent closure schemes (\forcode{ln_zdf{tke,gls,osm}})}
\label{subsec:ZDF_tcs}

The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and
\autoref{subsec:ZDF_osm} (\ie\ \np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} or \np{ln_zdfosm}{ln\_zdfosm} defined) deal, in theory,
with statically unstable density profiles.
In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in
\autoref{eq:ZDF_tke_e} or \autoref{eq:ZDF_gls_e} becomes a source term, since $N^2$ is negative.
It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also of the four neighboring values at
velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$).
These large values restore the static stability of the water column in a way similar to that of
the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}).
However, in the vicinity of the sea surface (first ocean layer), the eddy coefficients computed by
the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$,
because the mixing length scale is bounded by the distance to the sea surface.
It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme,
\ie\ setting the \np{ln_zdfnpc}{ln\_zdfnpc} namelist parameter to true and
defining the turbulent closure (\np{ln_zdftke}{ln\_zdftke} or \np{ln_zdfgls}{ln\_zdfgls} = \forcode{.true.}) all together.

The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection,
%as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
therefore \np[=.false.]{ln_zdfevd}{ln\_zdfevd} should be used with the OSMOSIS scheme.
% gm%  + one word on non local flux with KPP scheme trakpp.F90 module...

%% =================================================================================================
\section[Double diffusion mixing (\forcode{ln_zdfddm})]{Double diffusion mixing (\protect\np{ln_zdfddm}{ln\_zdfddm})}
\label{subsec:ZDF_ddm}

%\nlst{namzdf_ddm}

This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter
\np{ln_zdfddm}{ln\_zdfddm} in \nam{zdf}{zdf}.
Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa.
The former condition leads to salt fingering and the latter to diffusive convection.
Double-diffusive phenomena contribute to diapycnal mixing in extensive regions of the ocean.
\citet{merryfield.holloway.ea_JPO99} include a parameterisation of such phenomena in a global ocean model and show that
it leads to relatively minor changes in circulation but exerts significant regional influences on
temperature and salinity.

Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
\begin{align*}
  % \label{eq:ZDF_ddm_Kz}
  &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT}	\\
  &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
\end{align*}
where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
and $o$ by processes other than double diffusion.
The rates of double-diffusive mixing depend on the buoyancy ratio
$R_\rho = \alpha \partial_z T / \beta \partial_z S$, where $\alpha$ and $\beta$ are coefficients of
thermal expansion and saline contraction (see \autoref{subsec:TRA_eos}).
To represent mixing of $S$ and $T$ by salt fingering, we adopt the diapycnal diffusivities suggested by Schmitt
(1981):
\begin{align}
  \label{eq:ZDF_ddm_f}
  A_f^{vS} &=
             \begin{cases}
               \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
               0 				  					    &\text{otherwise}
             \end{cases}
  \\ 		    \label{eq:ZDF_ddm_f_T}
  A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho
\end{align}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_ddm}
  \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and
  diffusive convection]{
    From \citet{merryfield.holloway.ea_JPO99}:
    (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in
    regions of salt fingering.
    Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and
    thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
    (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in
    regions of diffusive convection.
    Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation.
    The latter is not implemented in \NEMO.}
  \label{fig:ZDF_ddm}
\end{figure}

The factor 0.7 in \autoref{eq:ZDF_ddm_f_T} reflects the measured ratio $\alpha F_T /\beta F_S \approx  0.7$ of
buoyancy flux of heat to buoyancy flux of salt (\eg, \citet{mcdougall.taylor_JMR84}).
Following  \citet{merryfield.holloway.ea_JPO99}, we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.

To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by
Federov (1988) is used:
\begin{align}
  % \label{eq:ZDF_ddm_d}
  A_d^{vT} &=
             \begin{cases}
               1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
               &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
               0 								&\text{otherwise}
             \end{cases}
                                       \nonumber \\
  \label{eq:ZDF_ddm_d_S}
  A_d^{vS} &=
             \begin{cases}
               A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
               A_d^{vT} \ 0.15 \ R_\rho               &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
               0 								&\text{otherwise}
             \end{cases}
\end{align}

The dependencies of \autoref{eq:ZDF_ddm_f} to \autoref{eq:ZDF_ddm_d_S} on $R_\rho$ are illustrated in
\autoref{fig:ZDF_ddm}.
Implementing this requires computing $R_\rho$ at each grid point on every time step.
This is done in \mdl{eosbn2} at the same time as $N^2$ is computed.
This avoids duplication in the computation of $\alpha$ and $\beta$ (which is usually quite expensive).

%% =================================================================================================
\section[Bottom and top friction (\textit{zdfdrg.F90})]{Bottom and top friction (\protect\mdl{zdfdrg})}
\label{sec:ZDF_drg}

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

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

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

Options to define the top and bottom friction are defined through the \nam{drg}{drg} namelist variables.
The bottom friction represents the friction generated by the bathymetry.
The top friction represents the friction generated by the ice shelf/ocean interface.
As the friction processes at the top and the bottom are treated in and identical way,
the description below considers mostly the bottom friction case, if not stated otherwise.

Both the surface momentum flux (wind stress) and the bottom momentum flux (bottom friction) enter the equations as
a condition on the vertical diffusive flux.
For the bottom boundary layer, one has:
 \[
   % \label{eq:ZDF_bfr_flux}
   A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
 \]
where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside
the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean).
How ${\cal F}_h^{\textbf U}$ influences the interior depends on the vertical resolution of the model near
the bottom relative to the Ekman layer depth.
For example, in order to obtain an Ekman layer depth $d = \sqrt{2\;A^{vm}} / f = 50$~m,
one needs a vertical diffusion coefficient $A^{vm} = 0.125$~m$^2$s$^{-1}$
(for a Coriolis frequency $f = 10^{-4}$~m$^2$s$^{-1}$).
With a background diffusion coefficient $A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.
When the vertical mixing coefficient is this small, using a flux condition is equivalent to
entering the viscous forces (either wind stress or bottom friction) as a body force over the depth of the top or
bottom model layer.
To illustrate this, consider the equation for $u$ at $k$, the last ocean level:
\begin{equation}
  \label{eq:ZDF_drg_flux2}
  \frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}}
\end{equation}
If the bottom layer thickness is 200~m, the Ekman transport will be distributed over that depth.
On the other hand, if the vertical resolution is high (1~m or less) and a turbulent closure model is used,
the turbulent Ekman layer will be represented explicitly by the model.
However, the logarithmic layer is never represented in current primitive equation model applications:
it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $.
Two choices are available in \NEMO: a linear and a quadratic bottom friction.
Note that in both cases, the rotation between the interior velocity and the bottom friction is neglected in
the present release of \NEMO.

In the code, the bottom friction is imposed by adding the trend due to the bottom friction to
 the general momentum trend in \mdl{dynzdf}.
For the time-split surface pressure gradient algorithm, the momentum trend due to
the barotropic component needs to be handled separately.
For this purpose it is convenient to compute and store coefficients which can be simply combined with
bottom velocities and geometric values to provide the momentum trend due to bottom friction.
 These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where:
\begin{equation}
  \label{eq:ZDF_bfr_bdef}
  \frac{\partial {\textbf U_h}}{\partial t} =
  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
\end{equation}
where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
Note than from \NEMO\ 4.0, drag coefficients are only computed at cell centers (\ie\ at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points.

%% =================================================================================================
\subsection[Linear top/bottom friction (\forcode{ln_lin})]{Linear top/bottom friction (\protect\np{ln_lin}{ln\_lin})}
\label{subsec:ZDF_drg_linear}

The linear friction parameterisation (including the special case of a free-slip condition) assumes that
the friction is proportional to the interior velocity (\ie\ the velocity of the first/last model level):
\[
  % \label{eq:ZDF_bfr_linear}
  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
\]
where $r$ is a friction coefficient expressed in $m s^{-1}$.
This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
and setting $r = H / \tau$, where $H$ is the ocean depth.
Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}.
A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used in quasi-geostrophic models.
One may consider the linear friction as an approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$
(\citet{gill_bk82}, Eq. 9.6.6).
For example, with a drag coefficient $C_D = 0.002$, a typical speed of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$,
and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days.
It can be changed by specifying \np{rn_Uc0}{rn\_Uc0} (namelist parameter).

 For the linear friction case the drag coefficient used in the general expression \autoref{eq:ZDF_bfr_bdef} is:
\[
  % \label{eq:ZDF_bfr_linbfr_b}
    c_b^T = - r
\]
When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}.
Setting \np[=.true.]{ln_drg_OFF}{ln\_drg\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.

These values are assigned in \mdl{zdfdrg}.
Note that there is support for local enhancement of these values via an externally defined 2D mask array
(\np[=.true.]{ln_boost}{ln\_boost}) given in the \textit{bfr\_coef.nc} input NetCDF file.
The mask values should vary from 0 to 1.
Locations with a non-zero mask value will have the friction coefficient increased by
$mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.

%% =================================================================================================
\subsection[Non-linear top/bottom friction (\forcode{ln_non_lin})]{Non-linear top/bottom friction (\protect\np{ln_non_lin}{ln\_non\_lin})}
\label{subsec:ZDF_drg_nonlinear}

The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic:
\[
  % \label{eq:ZDF_drg_nonlinear}
  {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h
  }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b
\]
where $C_D$ is a drag coefficient, and $e_b $ a top/bottom turbulent kinetic energy due to tides,
internal waves breaking and other short time scale currents.
A typical value of the drag coefficient is $C_D = 10^{-3} $.
As an example, the CME experiment \citep{treguier_JGR92} uses $C_D = 10^{-3}$ and
$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and
$e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
The CME choices have been set as default values (\np{rn_Cd0}{rn\_Cd0} and \np{rn_ke0}{rn\_ke0} namelist parameters).

As for the linear case, the friction is imposed in the code by adding the trend due to
the friction to the general momentum trend in \mdl{dynzdf}.
For the non-linear friction case the term computed in \mdl{zdfdrg} is:
\[
  % \label{eq:ZDF_drg_nonlinbfr}
    c_b^T = - \; C_D\;\left[ \left(\bar{u_b}^{i}\right)^2 + \left(\bar{v_b}^{j}\right)^2 + e_b \right]^{1/2}
\]

The coefficients that control the strength of the non-linear friction are initialised as namelist parameters:
$C_D$= \np{rn_Cd0}{rn\_Cd0}, and $e_b$ =\np{rn_bfeb2}{rn\_bfeb2}.
Note that for applications which consider tides explicitly, a low or even zero value of \np{rn_bfeb2}{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array
(\np[=.true.]{ln_boost}{ln\_boost}).
This works in the same way as for the linear friction case with non-zero masked locations increased by
$mask\_value$ * \np{rn_boost}{rn\_boost} * \np{rn_Cd0}{rn\_Cd0}.

%% =================================================================================================
\subsection[Log-layer top/bottom friction (\forcode{ln_loglayer})]{Log-layer top/bottom friction (\protect\np{ln_loglayer}{ln\_loglayer})}
\label{subsec:ZDF_drg_loglayer}

In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using
a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so.
If  \np[=.true.]{ln_loglayer}{ln\_loglayer}, $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):
\[
  C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2
\]

\noindent where $\kappa$ is the von-Karman constant and \np{rn_z0}{rn\_z0} is a roughness length provided via the namelist.

The drag coefficient is bounded such that it is kept greater or equal to
the base \np{rn_Cd0}{rn\_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter:
\np{rn_Cdmax}{rn\_Cdmax}, \ie
\[
  rn\_Cd0 \leq C_D \leq rn\_Cdmax
\]

\noindent The log-layer enhancement can also be applied to the top boundary friction if
under ice-shelf cavities are activated (\np[=.true.]{ln_isfcav}{ln\_isfcav}).
%In this case, the relevant namelist parameters are \np{rn_tfrz0}{rn\_tfrz0}, \np{rn_tfri2}{rn\_tfri2} and \np{rn_tfri2_max}{rn\_tfri2\_max}.

%% =================================================================================================
\subsection[Explicit top/bottom friction (\forcode{ln_drgimp=.false.})]{Explicit top/bottom friction (\protect\np[=.false.]{ln_drgimp}{ln\_drgimp})}
\label{subsec:ZDF_drg_stability}

Setting \np[=.false.]{ln_drgimp}{ln\_drgimp} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads:

At the top (below an ice shelf cavity):
\[
  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t}
  = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t}
\]

At the bottom (above the sea floor):
\[
  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b}
  = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b}
\]

Since this is conditionally stable, some care needs to exercised over the choice of parameters to ensure that the implementation of explicit top/bottom friction does not induce numerical instability.
For the purposes of stability analysis, an approximation to \autoref{eq:ZDF_drg_flux2} is:
\begin{equation}
  \label{eq:ZDF_Eqn_drgstab}
  \begin{split}
    \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
    &= -\frac{ru}{e_{3u}}\;2\rdt\\
  \end{split}
\end{equation}
\noindent where linear friction and a leapfrog timestep have been assumed.
To ensure that the friction cannot reverse the direction of flow it is necessary to have:
\[
  |\Delta u| < \;|u|
\]
\noindent which, using \autoref{eq:ZDF_Eqn_drgstab}, gives:
\[
  r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
\]
This same inequality can also be derived in the non-linear bottom friction case if
a velocity of 1 m.s$^{-1}$ is assumed.
Alternatively, this criterion can be rearranged to suggest a minimum bottom box thickness to ensure stability:
\[
  e_{3u} > 2\;r\;\rdt
\]
\noindent which it may be necessary to impose if partial steps are being used.
For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m.
For most applications, with physically sensible parameters these restrictions should not be of concern.
But caution may be necessary if attempts are made to locally enhance the bottom friction parameters.
To ensure stability limits are imposed on the top/bottom friction coefficients both
during initialisation and at each time step.
Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
The number of breaches of the stability criterion are reported as well as
the minimum and maximum values that have been set.
The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}.
Values of the friction coefficient are reduced as necessary to ensure stability;
these changes are not reported.

Limits on the top/bottom friction coefficient are not imposed if the user has elected to
handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}).
The number of potential breaches of the explicit stability criterion are still reported for information purposes.

%% =================================================================================================
\subsection[Implicit top/bottom friction (\forcode{ln_drgimp=.true.})]{Implicit top/bottom friction (\protect\np[=.true.]{ln_drgimp}{ln\_drgimp})}
\label{subsec:ZDF_drg_imp}

An optional implicit form of bottom friction has been implemented to improve model stability.
We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting.
This option can be invoked by setting \np{ln_drgimp}{ln\_drgimp} to \forcode{.true.} in the \nam{drg}{drg} namelist.
%This option requires \np{ln_zdfexp}{ln\_zdfexp} to be \forcode{.false.} in the \nam{zdf}{zdf} namelist.

This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step:

At the top (below an ice shelf cavity):
\[
  % \label{eq:ZDF_dynZDF__drg_top}
  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t}
  = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t}
\]

At the bottom (above the sea floor):
\[
  % \label{eq:ZDF_dynZDF__drg_bot}
  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b}
  = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b}
\]

where $t$ and $b$ refers to top and bottom layers respectively.
Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit.

%% =================================================================================================
\subsection[Bottom friction with split-explicit free surface]{Bottom friction with split-explicit free surface}
\label{subsec:ZDF_drg_ts}

With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \np[=.false.]{ln_drgimp}{ln\_drgimp} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie\ \forcode{ln_bt_fw=.false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \np[=.true.]{ln_drgimp}{ln\_drgimp},  stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.

The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO\ is as follows:
\begin{enumerate}
\item To extend the stability of the barotropic sub-stepping, bottom stresses are refreshed at each sub-iteration. The baroclinic part of the flow entering the stresses is frozen at the initial time of the barotropic iteration. In case of non-linear friction, the drag coefficient is also constant.
\item In case of an implicit drag, specific computations are performed in \mdl{dynzdf} which renders the overall scheme mixed explicit/implicit: the barotropic components of 3d velocities are removed before seeking for the implicit vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for thanks to the updated values of barotropic velocities. Then the implicit solution of 3d velocities is obtained. Lastly, the residual barotropic component is replaced by the time split estimate.
\end{enumerate}

Note that other strategies are possible, like considering vertical diffusion step in advance, \ie\ prior barotropic integration.

%% =================================================================================================
\section[Internal wave-driven mixing (\forcode{ln_zdfiwm})]{Internal wave-driven mixing (\protect\np{ln_zdfiwm}{ln\_zdfiwm})}
\label{subsec:ZDF_tmx_new}

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

The parameterization of mixing induced by breaking internal waves is a generalization of
the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}.
A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
and the resulting diffusivity is obtained as
\[
  % \label{eq:ZDF_Kwave}
  A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
\]
where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of
the energy available for mixing.
If the \np{ln_mevar}{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and
equal to 1/6 \citep{osborn_JPO80}.
In the opposite (recommended) case, $R_f$ is instead a function of
the turbulence intensity parameter $Re_b = \frac{ \epsilon}{\nu \, N^2}$,
with $\nu$ the molecular viscosity of seawater, following the model of \cite{bouffard.boegman_DAO13} and
the implementation of \cite{de-lavergne.madec.ea_JPO16}.
Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when
the mixing efficiency is constant.

In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
as a function of $Re_b$ by setting the \np{ln_tsdiff}{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice.
This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14},
is implemented as in \cite{de-lavergne.madec.ea_JPO16}.

The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$,
is constructed from three static maps of column-integrated internal wave energy dissipation,
$E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures:

\begin{align*}
  F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
  F_{pyc}(i,j,k) &\propto N^{n_p}\\
  F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
\end{align*}
In the above formula, $h_{ab}$ denotes the height above bottom,
$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
\[
  h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
\]
The $n_p$ parameter (given by \np{nn_zpyc}{nn\_zpyc} in \nam{zdf_iwm}{zdf\_iwm} namelist)
controls the stratification-dependence of the pycnocline-intensified dissipation.
It can take values of $1$ (recommended) or $2$.
Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of
the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.
$h_{cri}$ is related to the large-scale topography of the ocean (etopo2) and
$h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
the abyssal hill topography \citep{goff_JGR10} and the latitude.
% Jc: input files names ?

%% =================================================================================================
\section[Surface wave-induced mixing (\forcode{ln_zdfswm})]{Surface wave-induced mixing (\protect\np{ln_zdfswm}{ln\_zdfswm})}
\label{subsec:ZDF_swm}

Surface waves produce an enhanced mixing through wave-turbulence interaction.
In addition to breaking waves induced turbulence (\autoref{subsec:ZDF_tke}),
the influence of non-breaking waves can be accounted introducing
wave-induced viscosity and diffusivity as a function of the wave number spectrum.
Following \citet{qiao.yuan.ea_OD10}, a formulation of wave-induced mixing coefficient
is provided  as a function of wave amplitude, Stokes Drift and wave-number:

\begin{equation}
  \label{eq:ZDF_Bv}
  B_{v} = \alpha {A} {U}_{st} {exp(3kz)}
\end{equation}

Where $B_{v}$ is the wave-induced mixing coefficient, $A$ is the wave amplitude,
${U}_{st}$ is the Stokes Drift velocity, $k$ is the wave number and $\alpha$
is a constant which should be determined by observations or
numerical experiments and is set to be 1.

The coefficient $B_{v}$ is then directly added to the vertical viscosity
and diffusivity coefficients.

In order to account for this contribution set: \forcode{ln_zdfswm=.true.},
then wave interaction has to be activated through \forcode{ln_wave=.true.},
the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.}
(see \autoref{subsec:SBC_wave_sdw})
and the needed wave fields (significant wave height and mean wave number) can be provided either in forcing or coupled mode
(for more information on wave parameters and settings see \autoref{sec:SBC_wave})

%% =================================================================================================
\section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp})]{Adaptive-implicit vertical advection(\protect\np{ln_zad_Aimp}{ln\_zad\_Aimp})}
\label{subsec:ZDF_aimp}

The adaptive-implicit vertical advection option in NEMO is based on the work of
\citep{shchepetkin_OM15}.  In common with most ocean models, the timestep used with NEMO
needs to satisfy multiple criteria associated with different physical processes in order
to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical
CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the
constraints for a range of time and space discretizations and provide the CFL stability
criteria for a range of advection schemes. The values for the Leap-Frog with Robert
asselin filter time-stepping (as used in NEMO) are reproduced in
\autoref{tab:ZDF_zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these
restrictions but at the cost of large dispersive errors and, possibly, large numerical
viscosity. The adaptive-implicit vertical advection option provides a targetted use of the
implicit scheme only when and where potential breaches of the vertical CFL condition
occur. In many practical applications these events may occur remote from the main area of
interest or due to short-lived conditions such that the extra numerical diffusion or
viscosity does not greatly affect the overall solution. With such applications, setting:
\forcode{ln_zad_Aimp=.true.} should allow much longer model timesteps to be used whilst
retaining the accuracy of the high order explicit schemes over most of the domain.

\begin{table}[htbp]
  \centering
  % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}}
  \begin{tabular}{r|ccc}
    \hline
    spatial discretization  & 2$^nd$ order centered & 3$^rd$ order upwind & 4$^th$ order compact \\
    advective CFL criterion &                 0.904 &              0.472  &                0.522 \\
    \hline
  \end{tabular}
  \caption[Advective CFL criteria for the leapfrog with Robert Asselin filter time-stepping]{
    The advective CFL criteria for a range of spatial discretizations for
    the leapfrog with Robert Asselin filter time-stepping
    ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}.}
  \label{tab:ZDF_zad_Aimp_CFLcrit}
\end{table}

In particular, the advection scheme remains explicit everywhere except where and when
local vertical velocities exceed a threshold set just below the explicit stability limit.
Once the threshold is reached a tapered transition towards an implicit scheme is used by
partitioning the vertical velocity into a part that can be treated explicitly and any
excess that must be treated implicitly. The partitioning is achieved via a Courant-number
dependent weighting algorithm as described in \citep{shchepetkin_OM15}.

The local cell Courant number ($Cu$) used for this partitioning is:

\begin{equation}
  \label{eq:ZDF_Eqn_zad_Aimp_Courant}
  \begin{split}
    Cu &= {2 \rdt \over e^n_{3t_{ijk}}} \bigg (\big [ \texttt{Max}(w^n_{ijk},0.0) - \texttt{Min}(w^n_{ijk+1},0.0) \big ]    \\
       &\phantom{=} +\big [ \texttt{Max}(e_{{2_u}ij}e^n_{{3_{u}}ijk}u^n_{ijk},0.0) - \texttt{Min}(e_{{2_u}i-1j}e^n_{{3_{u}}i-1jk}u^n_{i-1jk},0.0) \big ]
                     \big / e_{{1_t}ij}e_{{2_t}ij}            \\
       &\phantom{=} +\big [ \texttt{Max}(e_{{1_v}ij}e^n_{{3_{v}}ijk}v^n_{ijk},0.0) - \texttt{Min}(e_{{1_v}ij-1}e^n_{{3_{v}}ij-1k}v^n_{ij-1k},0.0) \big ]
                     \big / e_{{1_t}ij}e_{{2_t}ij} \bigg )    \\
  \end{split}
\end{equation}

\noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as:

\begin{align}
  \label{eq:ZDF_Eqn_zad_Aimp_partition}
Cu_{min} &= 0.15 \nonumber \\
Cu_{max} &= 0.3  \nonumber \\
Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\
Fcu    &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\
\cf &=
     \begin{cases}
        0.0                                                        &\text{if $Cu \leq Cu_{min}$} \\
        (Cu - Cu_{min})^2 / (Fcu +  (Cu - Cu_{min})^2)             &\text{else if $Cu < Cu_{cut}$} \\
        (Cu - Cu_{max}) / Cu                                       &\text{else}
     \end{cases}
\end{align}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_coeff}
  \caption[Partitioning coefficient used to partition vertical velocities into parts]{
    The value of the partitioning coefficient (\cf) used to partition vertical velocities into
    parts to be treated implicitly and explicitly for a range of typical Courant numbers
    (\forcode{ln_zad_Aimp=.true.}).}
  \label{fig:ZDF_zad_Aimp_coeff}
\end{figure}

\noindent The partitioning coefficient is used to determine the part of the vertical
velocity that must be handled implicitly ($w_i$) and to subtract this from the total
vertical velocity ($w_n$) to leave that which can continue to be handled explicitly:

\begin{align}
  \label{eq:ZDF_Eqn_zad_Aimp_partition2}
    w_{i_{ijk}} &= \cf_{ijk} w_{n_{ijk}}     \nonumber \\
    w_{n_{ijk}} &= (1-\cf_{ijk}) w_{n_{ijk}}
\end{align}

\noindent Note that the coefficient is such that the treatment is never fully implicit;
the three cases from \autoref{eq:ZDF_Eqn_zad_Aimp_partition} can be considered as:
fully-explicit; mixed explicit/implicit and mostly-implicit.  With the settings shown the
coefficient (\cf) varies as shown in \autoref{fig:ZDF_zad_Aimp_coeff}. Note with these values
the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly
implicit' is 0.45 which is just below the stability limited given in
\autoref{tab:ZDF_zad_Aimp_CFLcrit}  for a 3rd order scheme.

The $w_i$ component is added to the implicit solvers for the vertical mixing in
\mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}.  This is
sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further
intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}).
For these schemes the implicit upstream fluxes must be added to both the monotonic guess
and to the higher order solution when calculating the antidiffusive fluxes. The implicit
vertical fluxes are then removed since they are added by the implicit solver later on.

The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be
used in a wide range of simulations. The following test simulation, however, does illustrate
the potential benefits and will hopefully encourage further testing and feedback from users:

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_frames}
  \caption[OVERFLOW: time-series of temperature vertical cross-sections]{
    A time-series of temperature vertical cross-sections for the OVERFLOW test case.
    These results are for the default settings with \forcode{nn_rdt=10.0} and
    without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}).}
  \label{fig:ZDF_zad_Aimp_overflow_frames}
\end{figure}

%% =================================================================================================
\subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case}

The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case}
provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case
by only a few extra physics choices namely:

\begin{forlines}
     ln_dynldf_OFF = .false.
     ln_dynldf_lap = .true.
     ln_dynldf_hor = .true.
     ln_zdfnpc     = .true.
     ln_traadv_fct = .true.
        nn_fct_h   =  2
        nn_fct_v   =  2
\end{forlines}

\noindent which were chosen to provide a slightly more stable and less noisy solution. The
result when using the default value of \forcode{nn_rdt=10.} without adaptive-implicit
vertical velocity is illustrated in \autoref{fig:ZDF_zad_Aimp_overflow_frames}. The mass of
cold water, initially sitting on the shelf, moves down the slope and forms a
bottom-trapped, dense plume. Even with these extra physics choices the model is close to
stability limits and attempts with \forcode{nn_rdt=30.} will fail after about 5.5 hours
with excessively high horizontal velocities. This time-scale corresponds with the time the
plume reaches the steepest part of the topography and, although detected as a horizontal
CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good
candidate, therefore, for use of the adaptive-implicit vertical advection scheme.

The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps
are shown in \autoref{fig:ZDF_zad_Aimp_overflow_all_rdt} (together with the equivalent
frames from the base run).  In this simple example the use of the adaptive-implicit
vertcal advection scheme has enabled a 12x increase in the model timestep without
significantly altering the solution (although at this extreme the plume is more diffuse
and has not travelled so far).  Notably, the solution with and without the scheme is
slightly different even with \forcode{nn_rdt=10.}; suggesting that the base run was
close enough to instability to trigger the scheme despite completing successfully.
To assist in diagnosing how active the scheme is, in both location and time, the 3D
implicit and explicit components of the vertical velocity are available via XIOS as
\texttt{wimp} and \texttt{wexp} respectively.  Likewise, the partitioning coefficient
(\cf) is also available as \texttt{wi\_cff}. For a quick oversight of
the schemes activity the global maximum values of the absolute implicit component
of the vertical velocity and the partitioning coefficient are written to the netCDF
version of the run statistics file (\texttt{run.stat.nc}) if this is active (see
\autoref{sec:MISC_opt} for activation details).

\autoref{fig:ZDF_zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for
the various overflow tests.  Note that the adaptive-implicit vertical advection scheme is
active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the
test case is close to stability limits even with this value. At the larger timesteps, the
vertical velocity is treated mostly implicitly at some location throughout the run. The
oscillatory nature of this measure appears to be linked to the progress of the plume front
as each cusp is associated with the location of the maximum shifting to the adjacent cell.
This is illustrated in \autoref{fig:ZDF_zad_Aimp_maxCf_loc} where the i- and k- locations of the
maximum have been overlaid for the base run case.

\medskip
\noindent Only limited tests have been performed in more realistic configurations. In the
ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes
restartability and reproducibility tests but it is unable to improve the model's stability
enough to allow an increase in the model time-step. A view of the time-series of maximum
partitioning coefficient (not shown here)  suggests that the default time-step of 5400s is
already pushing at stability limits, especially in the initial start-up phase. The
time-series does not, however, exhibit any of the 'cuspiness' found with the overflow
tests.

\medskip
\noindent A short test with an eORCA1 configuration promises more since a test using a
time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the
time-step is limited to 2700s without.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_all_rdt}
  \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{
    Sample temperature vertical cross-sections from mid- and end-run using
    different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection.
    Without the adaptive implicit vertical advection
    only the run with the shortest timestep is able to run to completion.
    Note also that the colour-scale has been chosen to confirm that
    temperatures remain within the original range of 10$^o$ to 20$^o$.}
  \label{fig:ZDF_zad_Aimp_overflow_all_rdt}
\end{figure}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf}
  \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{
    The maximum partitioning coefficient during a series of test runs with
    increasing model timestep length.
    At the larger timesteps,
    the vertical velocity is treated mostly implicitly at some location throughout the run.}
  \label{fig:ZDF_zad_Aimp_maxCf}
\end{figure}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf_loc}
  \caption[OVERFLOW: maximum partitioning coefficient for the case overlaid]{
    The maximum partitioning coefficient for the \forcode{nn_rdt=10.0} case overlaid with
    information on the gridcell i- and k-locations of the maximum value.}
  \label{fig:ZDF_zad_Aimp_maxCf_loc}
\end{figure}

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

\end{document}