chap_ZDF.tex

vertical diffusion (\autoref{subsec:ZDF_evd}).
It is recommended that one of these parameterisations be enabled when using either the TKE or GLS turbulent
closure scheme, but not when using the OSMOSIS turbulent closure scheme (see \autoref{subsec:ZDF_tcs}).

%% =================================================================================================
\subsection[Non-penetrative convective adjustment (\forcode{ln_tranpc})]{Non-penetrative convective adjustment (\protect\np{ln_tranpc}{ln\_tranpc})}
\label{subsec:ZDF_npc}

\begin{figure}[!htb]
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_npc}
  \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{
    Example of an unstable density profile treated by
    the non penetrative convective adjustment algorithm.
    $1^{st}$ step: the initial profile is checked from the surface to the bottom.
    It is found to be unstable between levels 3 and 4.
    They are mixed.
    The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed.
    The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed.
    The $1^{st}$ step ends since the density profile is then stable below the level 3.
    $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step:
    levels 2 to 5 are mixed.
    The new density profile is checked.
    It is found stable: end of algorithm.}
  \label{fig:ZDF_npc}
\end{figure}

Options are defined through the \nam{zdf}{zdf} namelist variables.
The non-penetrative convective adjustment is used when \np[=.true.]{ln_zdfnpc}{ln\_zdfnpc}.
It is applied at each \np{nn_npc}{nn\_npc} time step and mixes downwards instantaneously the statically unstable portion of
the water column, but only until the density structure becomes neutrally stable
(\ie\ until the mixed portion of the water column has \textit{exactly} the density of the water just below)
\citep{madec.delecluse.ea_JPO91}.
The associated algorithm is an iterative process used in the following way (\autoref{fig:ZDF_npc}):
starting from the top of the ocean, the first instability is found.
Assume in the following that the instability is located between levels $k$ and $k+1$.
The temperature and salinity in the two levels are vertically mixed, conserving the heat and salt contents of
the water column.
The new density is then computed by a linear approximation.
If the new density profile is still unstable between levels $k+1$ and $k+2$,
levels $k$, $k+1$ and $k+2$ are then mixed.
This process is repeated until stability is established below the level $k$
(the mixing process can go down to the ocean bottom).
The algorithm is repeated to check if the density profile between level $k-1$ and $k$ is unstable and/or
if there is no deeper instability.

This algorithm is significantly different from mixing statically unstable levels two by two.
The latter procedure cannot converge with a finite number of iterations for some vertical profiles while
the algorithm used in \NEMO\ converges for any profile in a number of iterations which is less than
the number of vertical levels.
This property is of paramount importance as pointed out by \citet{killworth_iprc89}:
it avoids the existence of permanent and unrealistic static instabilities at the sea surface.
This non-penetrative convective algorithm has been proved successful in studies of the deep water formation in
the north-western Mediterranean Sea \citep{madec.delecluse.ea_JPO91, madec.chartier.ea_DAO91, madec.crepon_iprc91}.

The current implementation has been modified in order to deal with any non linear equation of seawater
(L. Brodeau, personnal communication).
Two main differences have been introduced compared to the original algorithm:
$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency
(not the difference in potential density);
$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in
the same way their temperature and salinity has been mixed.
These two modifications allow the algorithm to perform properly and accurately with TEOS10 or EOS-80 without
having to recompute the expansion coefficients at each mixing iteration.

%% =================================================================================================
\subsection[Enhanced vertical diffusion (\forcode{ln_zdfevd})]{Enhanced vertical diffusion (\protect\np{ln_zdfevd}{ln\_zdfevd})}
\label{subsec:ZDF_evd}

Options are defined through the  \nam{zdf}{zdf} namelist variables.
The enhanced vertical diffusion parameterisation is used when \np[=.true.]{ln_zdfevd}{ln\_zdfevd}.
In this case, the vertical eddy mixing coefficients are assigned very large values
in regions where the stratification is unstable
(\ie\ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.
This is done either on tracers only (\np[=0]{nn_evdm}{nn\_evdm}) or
on both momentum and tracers (\np[=1]{nn_evdm}{nn\_evdm}).

In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and if \np[=1]{nn_evdm}{nn\_evdm},
the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ values also, are set equal to
the namelist parameter \np{rn_evd}{rn\_evd}.
A typical value for \np{rn_evd}{rn\_evd} is between 1 and $100~m^{2}s^{-1}$.
This parameterisation of convective processes is less time consuming than
the convective adjustment algorithm presented above when mixing both tracers and
momentum in the case of static instabilities.

Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$.
This removes a potential source of divergence of odd and even time step in
a leapfrog environment \citep{leclair_phd10} (see \autoref{sec:TD_mLF}).

%% =================================================================================================
\subsection[Mass Flux Convection (\forcode{ln_zdfmfc})]{Mass Flux Convection (\protect\np{ln_zdfmfc}{ln\_zdfmfc})}
\label{subsec:ZDF_mfc}

The \np{ln_zdfmfc}{ln\_zdfmfc} option offers a new, coherent way to simultaneously parameterize 
local and non-local transport within the oceanic convective mixing layer \citep{giordani.ea_JAMES2020}. 
This approach, initially developed for atmospheric models (\cite{grant_QJRMS01}; \cite{soares.ea_QJRMS04}; \cite{pergaud.ea_BLM09}), 
assumes that subgrid turbulent fluxes in the convective boundary layer result from two distinct mixing scales: 
the local scale of small eddies, represented by vertical diffusion (TKE, GLS, etc.), 
and the scale of large eddies or convective thermals, represented by a non-local mass flux approach. 
The combination of both diffusive and convective schemes operating simultaneously is called as Eddy Diffusivity Mass Flux (EDMF).

The mass flux scheme is designed to represent all convection regimes, i.e., from moderate to strong, and from intermediate to deep convection. 
Its goal is to calculate the effects of a population of convective plumes occupying a fraction of an ocean model grid cell on its prognostic variables. 
This scheme thus parameterizes the subgrid-scale effects of very fine-scale convective plumes on the resolved variables at the model grid scale. 
This phenomenon is a one-dimensional vertical process. Since convection is a non-local phenomenon, 
this scheme complements mixing schemes such as the TKE and Richardson schemes. 
It serves as an alternative to the "Enhanced Vertical Diffusion" (\np{ln_zdfevd}{ln\_zdfevd}) or
"Non-Penetrative Convection" (\np{ln_zdfnpc}{ln\_zdfnpc}) schemes currently implemented in NEMO. 
The scheme is active at every model time step, alongside diffusion. 
It is broken down into several steps, corresponding to the evolution equations 
for tracers within convective plumes and the mass flux resulting from the work of
buoyancy forces in the vertical direction and the conservation of energy. 
The system of equations, along with its implementation in NEMO, is described in
\citep{giordani.ea_JAMES2020}, specifically in equations 12 and 13.

%% =================================================================================================
\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 TKE and GLS turbulent closure schemes presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls}
are, in theory, able to handle 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}~ms^{-1}$,
because the mixing length scale is bounded by the distance to the sea surface.

When using either the TKE (\np[=.true.]{ln_zdftke}{ln\_zdftke}) or GLS (\np[=.true.]{ln_zdfgls}{ln\_zdfgls})
turbulent closure scheme, it is therefore recommended to also enable the enhanced vertical diffusion
parameterisation (\np[=.true.]{ln_zdfevd}{ln\_zdfevd}).
This should not be done when using the OSMOSIS turbulent closure scheme (\np[=.true.]{ln_zdfosm}{ln\_zdfosm}), as this
already includes enhanced vertical diffusion in the case of convection (see \autoref{subsec:ZDF_osm}).
%as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
% 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
\citet{federov1988_ddm} are 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} (\autoref{sec:TRA_eosbn2}) at the same time as $N^2$ is computed,
avoiding 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 via parameters in the \nam{drg}{drg} namelist,
and, for top and bottom friction specifically, in the \nam{drg_top}{drg\_top} and \nam{drg_bot}{drg\_bot}
namelists respectively.
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 an 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}$ 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 value, $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} (\autoref{sec:DYN_zdf}).
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 that from \NEMO\ 4.0, drag coefficients are only computed at cell centers (\ie\ at T-points) and are
referred 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[Free-slip boundary conditions (\forcode{ln_drg_OFF})]{Free-slip boundary conditions (\protect\np{ln_drg_OFF}{ln\_drg\_OFF})}
\label{subsec:ZDF_drg_free}

When setting \np[=.true.]{ln_drg_OFF}{ln\_drg\_OFF} free-slip conditions are used at the top and bottom boundaries,
\ie\ the drag coefficient used in \autoref{eq:ZDF_bfr_bdef} is set to $c_b^T = 0$.

%% =================================================================================================
\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, \autoref{subsec:ZDF_drg_free})
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}$.
The drag coefficient used in the general expression (\autoref{eq:ZDF_bfr_bdef}) is therefore:

\[
  % \label{eq:ZDF_bfr_linbfr_b}
    c_b^T = - r
\]

$r$ is generally estimated as $H / \tau$, where $\tau$ is a typical decay time in the deep ocean and $H$ is the ocean depth.
Commonly accepted values of $\tau$ are of the order of 100 to 200 days \citep{weatherly_JMR84}.
A value of $\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}$
\citep[][Eq. 9.6.6]{gill_bk82}. \\

In \NEMO, linear friction is enabled by setting \np[=.true.]{ln_lin}{ln\_lin}.
$c_b^T$ is calculated in \mdl{zdfdrg} and the trend due to the friction is added to the general momentum trend
in \mdl{dynzdf}.
$r$ is calculated as \forcode{rn_Cd0 * rn_Uc0}, where \np{rn_Cd0}{rn\_Cd0} and \np{rn_Uc0}{rn\_Uc0} are
namelist parameters corresponding to the drag coefficient $C_D$ and velocity scale $U_{av}$ respectively.
Their default values (0.001 and $0.4$~m\;s$^{-1}$ respectively) result in a friction coefficient of
$r = 4\;10^{-4}$~m\;s$^{-1}$, corresponding to a decay time scale of 115~days when assuming an ocean
depth of $H = 4000$~m.

Local enhancements may be applied to the values of $c_b^T$ by setting \np[=.true.]{ln_boost}{ln\_boost}
and providing a 2D mask array (with values $0 \leq M_b \leq 1$) via a NetCDF file.
For bottom (top) friction, the array and file are named \forcode{bfr_coef} and \textit{bfr\_coef.nc}
(\forcode{tfr_coef} and \textit{tfr\_coef.nc}) respectively.
Locations with a non-zero mask value will have the friction coefficient increased by
$M_b * \mathrm{rn\_boost} * \mathrm{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.
The drag coefficient used in the general expression (\autoref{eq:ZDF_bfr_bdef}) is therefore:
\[
  % \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}
\]

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}$. \\

In \NEMO, non-linear friction is enabled by setting \np[=.true.]{ln_non_lin}{ln\_non\_lin}.
As for linear friction, $c_b^T$ is calculated in \mdl{zdfdrg} and the trend due to the friction is
added to the general momentum trend in \mdl{dynzdf}.
$C_D$ and $e_b$ correspond to the namelist parameters \np{rn_Cd0}{rn\_Cd0} and \np{rn_ke0}{rn\_ke0} respectively,
with their default values set to those of the CME experiment.
Note that for applications which consider tides explicitly, a low or even zero value of \np{rn_ke0}{rn\_ke0} is recommended.

As for linear friction, local enhancements may be applied to the values of $c_b^T$ by setting
\np[=.true.]{ln_boost}{ln\_boost} (see \autoref{subsec:ZDF_drg_linear}).

%% =================================================================================================
\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.
$C_D$ is then no longer constant, but instead 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} / z_0 \right ) } \right )^2
\]

\noindent where $\kappa$ is the von-Karman constant and $z_0$ is a roughness length.
This assumes that the model vertical resolution can capture the logarithmic layer,
which typically occurs for layers thinner than 1 m or so.

This special case of the non-linear friction is enabled by setting
\np[=.true.]{ln_loglayer}{ln\_loglayer} instead of \np[=.true.]{ln_non_lin}{ln\_non\_lin}.
$z_0$ corresponds to the namelist parameter \np{rn_z0}{rn\_z0}, while $C_D$ is bounded by
the namelist parameters \np{rn_Cd0}{rn\_Cd0} and \np{rn_Cdmax}{rn\_Cdmax} such that
$\mathrm{rn\_Cd0} \leq C_D \leq \mathrm{rn\_Cdmax}$. The lower bound of \np{rn_Cd0}{rn\_Cd0}
covers large layer thicknesses where logarithmic layers are presumably not resolved, while the upper bound
of \np{rn_Cdmax}{rn\_Cdmax} is applied for stability reasons.

\noindent The log-layer enhancement can also be applied to the top boundary friction if
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~\mathrm{ms}^{-1}$, $\rdt = 1800~\mathrm{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~\mathrm{ms}^{-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:D2D_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 it difficult to obtain 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 obtaining the implicit
      vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for
      due 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 \citep{de-lavergne.vic.ea_JAMES20} is a generalization of the approach originally proposed by \citet{st-laurent.simmons.ea_GRL02}.
A three-dimensional field of turbulence production by breaking internal waves $\epsilon(i,j,k)$ 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 \citep{osborn_JPO80}.
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 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$^{-1}$, a limit that is often reached in weak stratification when the mixing efficiency is constant. 
In addition to the mixing efficiency, the ratio of salt to heat diffusivities can be 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}.

This parameterization of internal wave-driven mixing is energetically constrained because vertical diffusivities derive from a known and globally constant power supply \citep{de-lavergne.vic.ea_JAMES20}. In the routine \mdl{zdfiwm}, the background diffusivity (\forcode{avtb}) is forced to a very small value (10$^{-10}$ m$^2$ s$^{-1}$), irrespective of namelist choices for background mixing. In addition, in the \mdl{zdftke} routine, when \np{ln_zdfiwm}{ln\_zdfiwm} is set to \forcode{.true.}, the background level of turbulent kinetic energy and the minimum mixing length are set to small values avoiding background diffusivities higher than molecular rates of diffusion. Hence, when the parameterization of internal wave-driven mixing is activated, vertical mixing in the deep ocean is truly energy constrained.

The three-dimensional distribution of turbulence production, $\epsilon(i,j,k)$ [W kg$^{-1}$],
is constructed from four static maps of column-integrated internal wave energy dissipation [W m$^{-2}$] (\autoref{fig:ZDF_iwm}):
\begin{itemize}
\item $E_{nsq}(i,j)$: dissipation scaling with the squared buoyancy frequency $N^2$;
\item $E_{sho}(i,j)$: dissipation due to shoaling (scales with $N$ in the vertical);
\item $E_{cri}(i,j)$: bottom-intensified dissipation at critical slopes;
\item $E_{bot}(i,j)$: bottom-intensified dissipation above abyssal hills.
\end{itemize}
These power fields are read in a NetCDF forcing file whose name is set in the namelist section \nam{zdf}{zdf}. This forcing file also contains two maps of decay scales [m] necessary for the vertical distribution of $E_{cri}$ and $E_{bot}$:
\begin{itemize}
\item $H_{cri}$(i,j): related to the height difference of the critical slope.
\item $H_{bot}$(i,j): related to the wavelength of abyssal hills and to the energy flux $E_{bot}$.
\end{itemize}
Each power field goes with a specific vertical structure. The local turbulence production $\epsilon(i,j,k)$ is thus obtained as the sum of four contributions,
$$ \epsilon = \epsilon_{nsq} + \epsilon_{sho} + \epsilon_{cri} + \epsilon_{bot} ,$$
where, within each $(i,j)$ water column,
\begin{itemize}
\item $\epsilon_{nsq} \propto N^2$
\item $\epsilon_{sho} \propto N$
\item $\epsilon_{cri} \propto \exp(-h_{ab}/H_{cri})$
\item $\epsilon_{bot} \propto 1/(1+h_{ab}/H_{bot})^2$
\end{itemize}
with $h_{ab}$ the height above bottom and $N$ the simulated buoyancy frequency. These vertical structures are interactive in the sense that they depend on the simulated stratification and sea surface height. Vertical distributions are thus computed every time step and are normalized so that
$$\int \rho\,\epsilon_{nsq}\,dz = E_{nsq}\,,\; \int \rho\,\epsilon_{sho}\,dz = E_{sho}\,,\; \int \rho\,\epsilon_{cri}\,dz = E_{cri}\,,\; \int \rho\,\epsilon_{bot}\,dz = E_{bot}\;.$$

\begin{figure}
  \centering
  \includegraphics[width=0.66\textwidth]{ZDF_iwm}
  \caption[Internal tide energy dissipation]{
    Four power maps (in log$_{10}$(W m$^{-2}$)), as estimated by \cite{de-lavergne.rathore.ea_ESS24}, that enter the parameterization of internal wave-driven mixing. The globally integrated power (in GW) within each field is indicated at the top left of each panel. The overall power amounts to 1.40 TW.}
  \label{fig:ZDF_iwm}
\end{figure}

As of now, the parameterization only accounts for mixing powered by internal tides, thought to be the dominant supply of small-scale turbulence in the ocean interior \citep{de-lavergne.vic.ea_JAMES20}. It is recommended to use the forcing fields available at https://doi.org/10.17882/103233. These fields are state-of-the-art estimates of internal tide energy dissipation, including subinertial tides \citep{de-lavergne.rathore.ea_ESS24}. Details on the implementation of this parameterization and impacts in long NEMO eORCA1 experiments can be found in \citep{de-lavergne.rathore.ea_ESS24}.


%% =================================================================================================
\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 for by 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.

This parameterisation is enabled by setting \np[=.true.]{ln_zdfswm}{ln\_zdfswm}.
Additionally, both wave interaction (\np[=.true.]{ln_wave}{ln\_wave}) and calculation of the Stokes Drift
(\np[=.true.]{ln_sdw}{ln\_sdw}) must be enabled.
The required wave fields (significant wave height and mean wave number) can be provided either in forced or coupled mode.
For more information on wave parameters and settings, see \autoref{sec:SBC_wave}.


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

\end{document}