source: tags/nemo_v3_2/nemo_v3_2/DOC/TexFiles/Chapters/Chap_ZDF.tex @ 1878

% ================================================================
% Chapter Ñ Vertical Ocean Physics (ZDF)
% ================================================================
\chapter{Vertical Ocean Physics (ZDF)}
\label{ZDF}
\minitoc

%gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
\gmcomment{Steven remark (not taken into account : problem here with turbulent vs turbulence.  I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation!
Gurvan :  I kept "turbulent closure etc "...}


% ================================================================
% Vertical Mixing
% ================================================================
\section{Vertical Mixing}
\label{ZDF_zdf}

The discrete form of the ocean subgrid scale physics has been presented in 
\S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries, 
the turbulent fluxes of momentum, heat and salt have to be defined. At the 
surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}), 
while at the bottom they are set to zero for heat and salt, unless a geothermal 
flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} 
defined, see \S\ref{TRA_bbc}), and specified through a bottom friction 
parameterisation for momentum (see \S\ref{ZDF_bfr}).

In this section we briefly discuss the various choices offered to compute 
the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , 
$A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- 
points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These 
coefficients can be assumed to be either constant, or a function of the local 
Richardson number, or computed from a turbulent closure model (either 
TKE or KPP formulation). The computation of these coefficients is initialized 
in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or 
\mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer 
diffusion, including the surface forcing, are computed and added to the 
general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively. 
These trends can be computed using either a forward time stepping scheme 
(namelist parameter \np{np\_zdfexp}=true) or a backward time stepping 
scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing 
coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}).

% -------------------------------------------------------------------------------------------------------------
%        Constant 
% -------------------------------------------------------------------------------------------------------------
\subsection{Constant (\key{zdfcst})}
\label{ZDF_cst}
%--------------------------------------------namzdf---------------------------------------------------------
\namdisplay{namzdf}
%--------------------------------------------------------------------------------------------------------------

When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients 
are set to constant values over the whole ocean. This is the crudest way to define 
the vertical ocean physics. It is recommended that this option is only used in 
process studies, not in basin scale simulations. Typical values used in this case are:
\begin{align*} 
A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}  \\
\\
A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
\end{align*}

These values are set through the \np{avm0} and \np{avt0} namelist parameters. 
In all cases, do not use values smaller that those associated with the molecular 
viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, 
$\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity.


% -------------------------------------------------------------------------------------------------------------
%        Richardson Number Dependent
% -------------------------------------------------------------------------------------------------------------
\subsection{Richardson Number Dependent (\key{zdfric})}
\label{ZDF_ric}

%--------------------------------------------namric---------------------------------------------------------
\namdisplay{namric}
%--------------------------------------------------------------------------------------------------------------

When \key{zdfric} is defined, a local Richardson number dependent formulation 
for the vertical momentum and tracer eddy coefficients is set. The vertical mixing 
coefficients are diagnosed from the large scale variables computed by the model. 
\textit{In situ} measurements have been used to link vertical turbulent activity to 
large scale ocean structures. The hypothesis of a mixing mainly maintained by the 
growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 
vertical eddy coefficients and the local Richardson number ($i.e.$ the 
ratio of stratification to vertical shear). Following \citet{PacPhil1981}, the following 
formulation has been implemented:
\begin{equation} \label{Eq_zdfric}
   \left\{      \begin{aligned}
         A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
         \\
         A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
   \end{aligned}    \right.
\end{equation}
where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson 
number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
$A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the 
constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 
is the maximum value that can be reached by the coefficient when $Ri\leq 0$, 
$a=5$ and $n=2$. The last three values can be modified by setting the 
\np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively.

% -------------------------------------------------------------------------------------------------------------
%        TKE Turbulent Closure Scheme 
% -------------------------------------------------------------------------------------------------------------
\subsection{TKE Turbulent Closure Scheme (\key{zdftke})}
\label{ZDF_tke}

%--------------------------------------------namtke---------------------------------------------------------
\namdisplay{namtke}
%--------------------------------------------------------------------------------------------------------------

The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 
turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent 
kinetic energy, and a closure assumption for the turbulent length scales. This 
turbulent closure model has been developed by \citet{Bougeault1989} in the 
atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 
embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since 
then, significant modifications have been introduced by \citet{Madec1998} in both 
the implementation and the formulation of the mixing length scale. The time 
evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical 
shear, its destruction through stratification, its vertical diffusion, and its dissipation 
of \citet{Kolmogorov1942} type:
\begin{equation} \label{Eq_zdftke_e}
\frac{\partial \bar{e}}{\partial t} = 
\frac{A^{vm}}{e_3 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
                     +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
-A^{vT}\,N^2
+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
            \;\frac{\partial \bar{e}}{\partial k}} \right]
- c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
\end{equation}
\begin{equation} \label{Eq_zdftke_kz}
   \begin{split}
         A^{vm} &= C_k\  l_k\  \sqrt {\bar{e}}     \\
         A^{vT} &= A^{vm} / P_{rt}
   \end{split}
\end{equation}
where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 
$P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and 
$C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth 
\citep{Gaspar1990}. They are set through namelist parameters \np{ediff} 
and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993}, 
be a function of the local Richardson number, $R_i $:
\begin{align*} \label{Eq_prt}
P_{rt} = \begin{cases}
                    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}  \\
                    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}  \\
                    \ \ 10 &      \text{if $\ 2 \leq R_i$} 
            \end{cases}
\end{align*}
Note that a horizontal Shapiro filter can optionally be applied to $R_i$. 
However it is an obsolescent option that is not recommended.  
The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter.

For computational efficiency, the original formulation of the turbulent length 
scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 
are proposed, the choice of which is controlled by the \np{nmxl} namelist 
parameter. The first two are based on the following first order approximation 
\citep{Blanke1993}:
\begin{equation} \label{Eq_tke_mxl0_1}
l_k = l_\epsilon = \sqrt {2 \bar e} / N
\end{equation}
which is valid in a stable stratified region with constant values of the brunt-
Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance 
to the surface or to the bottom (\np{nmxl}=0) or by the local vertical scale factor (\np{nmxl}=1). \citet{Blanke1993} notice that this simplification has two major 
drawbacks: it makes no sense for locally unstable stratification and the 
computation no longer uses all the information contained in the vertical density 
profile. To overcome these drawbacks, \citet{Madec1998} introduces the 
\np{nmxl}=2 or 3 cases, which add an extra assumption concerning the vertical 
gradient of the computed length scale. So, the length scales are first evaluated 
as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
\begin{equation} \label{Eq_tke_mxl_constraint}
\frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
\qquad \text{with }\  l =  l_k = l_\epsilon
\end{equation}
\eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length 
scale cannot be larger than the variations of depth. It provides a better 
approximation of the \citet{Gaspar1990} formulation while being much less 
time consuming. In particular, it allows the length scale to be limited not only 
by the distance to the surface or to the ocean bottom but also by the distance 
to a strongly stratified portion of the water column such as the thermocline 
(Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} 
constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 
the upward and downward length scales, and evaluate the dissipation and 
mixing turbulent length scales as (and note that here we use numerical 
indexing):
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \label{Fig_mixing_length}  \begin{center}
\includegraphics[width=1.00\textwidth]{./TexFiles/Figures/Fig_mixing_length.pdf}
\caption {Illustration of the mixing length computation. }
\end{center}  
\end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{equation} \label{Eq_tke_mxl2}
\begin{aligned}
 l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)}\ \ \ \;  \right)
    \quad &\text{ from $k=1$ to $jpk$ }\ \\
 l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3T}^{(k-1)}  \right)   
    \quad &\text{ from $k=jpk$ to $1$ }\ \\
\end{aligned}
\end{equation}
where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1}, 
$i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$.

In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same 
value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the 
\np{nmxl}=2 case, the dissipation and mixing length scales are give 
as in \citet{Gaspar1990}:
\begin{equation} \label{Eq_tke_mxl_gaspar}
\begin{aligned}
& l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }    \\
& l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right) 
\end{aligned}
\end{equation}

At the sea surface the value of $\bar{e}$ is prescribed from the wind 
stress field: $\bar{e}=ebb\;\left| \tau \right|$ ($ebb=60$ by default) 
with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist
parameters). Its value at the bottom of the ocean is assumed to be 
equal to the value of the level just above. The time integration of the 
$\bar{e}$ equation may formally lead to negative values because the 
numerical scheme does not ensure its positivity. To overcome this 
problem, a cut-off in the minimum value of $\bar{e}$ is used. Following 
\citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. 
This allows the subsequent formulations to match that of\citet{Gargett1984} 
for the diffusion in the thermocline and deep ocean :  $(A^{vT} = 10^{-3} / N)$. 
In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical 
instabilities associated with too weak vertical diffusion. They must be 
specified at least larger than the molecular values, and are set through 
\textit{avm0} and \textit{avt0} (\textbf{namelist} parameters).

% -------------------------------------------------------------------------------------------------------------
%        K Profile Parametrisation (KPP) 
% -------------------------------------------------------------------------------------------------------------
\subsection{K Profile Parametrisation (KPP) (\key{zdfkpp}) }
\label{ZDF_kpp}

%--------------------------------------------namkpp--------------------------------------------------------
\namdisplay{namkpp}
%--------------------------------------------------------------------------------------------------------------

The K-Profile Parametrization (KKP) developed by \cite{Large_al_RG94} has been 
implemented in \NEMO by J. Chanut (PhD reference to be added here!).

\colorbox{yellow}{Add a description of KPP here.}


% ================================================================
% Convection
% ================================================================
\section{Convection}
\label{ZDF_conv}

%--------------------------------------------namzdf--------------------------------------------------------
\namdisplay{namzdf}
%--------------------------------------------------------------------------------------------------------------

Static instabilities (i.e. light potential densities under heavy ones) may 
occur at particular ocean grid points. In nature, convective processes 
quickly re-establish the static stability of the water column. These 
processes have been removed from the model via the hydrostatic 
assumption so they must be parameterized. Three parameterisations 
are available to deal with convective processes: a non-penetrative 
convective adjustment or an enhanced vertical diffusion, or/and the 
use of a turbulent closure scheme.

% -------------------------------------------------------------------------------------------------------------
%       Non-Penetrative Convective Adjustment 
% -------------------------------------------------------------------------------------------------------------
\subsection   [Non-Penetrative Convective Adjustment (\np{ln\_tranpc}) ]
         {Non-Penetrative Convective Adjustment (\np{ln\_tranpc}=.true.) }
\label{ZDF_npc}

%--------------------------------------------namnpc--------------------------------------------------------
\namdisplay{namnpc}
%--------------------------------------------------------------------------------------------------------------


%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!htb] \label{Fig_npc}   \begin{center}
\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_npc.pdf}
\caption {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.}
\end{center}   \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true. 
It is applied at each \np{nnpc1} time step and mixes downwards instantaneously 
the statically unstable portion of the water column, but only until the density 
structure becomes neutrally stable ($i.e.$ until the mixed portion of the water 
column has \textit{exactly} the density of the water just below) \citep{Madec1991a}. 
The associated algorithm is an iterative process used in the following way 
(Fig. \ref{Fig_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 potential 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{Killworth1989}: 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{Madec1991a, Madec1991b, Madec1991c}.

Note that in the current implementation of this algorithm presents several 
limitations. First, potential density referenced to the sea surface is used to 
check whether the density profile is stable or not. This is a strong 
simplification which leads to large errors for realistic ocean simulations. 
Indeed, many water masses of the world ocean, especially Antarctic Bottom
Water, are unstable when represented in surface-referenced potential density. 
The scheme will erroneously mix them up. Second, the mixing of potential 
density is assumed to be linear. This assures the convergence of the algorithm 
even when the equation of state is non-linear. Small static instabilities can thus 
persist due to cabbeling: they will be treated at the next time step. 
Third, temperature and salinity, and thus density, are mixed, but the 
corresponding velocity fields remain unchanged. When using a Richardson 
Number dependent eddy viscosity, the mixing of momentum is done through 
the vertical diffusion: after a static adjustment, the Richardson Number is zero 
and thus the eddy viscosity coefficient is at a maximum. When this convective 
adjustment algorithm is used with constant vertical eddy viscosity, spurious 
solutions can occur since the vertical momentum diffusion remains small even 
after a static adjustment. In that case, we recommend the addition of momentum 
mixing in a manner that mimics the mixing in temperature and salinity 
\citep{Speich1992, Speich1996}.

% -------------------------------------------------------------------------------------------------------------
%       Enhanced Vertical Diffusion 
% -------------------------------------------------------------------------------------------------------------
\subsection   [Enhanced Vertical Diffusion (\np{ln\_zdfevd})]
         {Enhanced Vertical Diffusion (\np{ln\_zdfevd}=.true.)}
\label{ZDF_evd}

%--------------------------------------------namzdf--------------------------------------------------------
\namdisplay{namzdf}
%--------------------------------------------------------------------------------------------------------------

The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true. 
In this case, the vertical eddy mixing coefficients are assigned very large values 
(a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 
($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}. 
This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and 
tracers (\np{n\_evdm}=1).

In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and 
if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
values also, are set equal to the namelist parameter \np{avevd}. A typical value 
for $avevd$ 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. It requires the use of an implicit time stepping on 
vertical diffusion terms (i.e. \np{ln\_zdfexp}=false). 

% -------------------------------------------------------------------------------------------------------------
%       Turbulent Closure Scheme 
% -------------------------------------------------------------------------------------------------------------
\subsection{Turbulent Closure Scheme (\key{zdftke})}
\label{ZDF_tcs}

The TKE turbulent closure scheme presented in \S\ref{ZDF_tke} and used 
when the \key{zdftke} is defined, in theory solves the problem of statically 
unstable density profiles. In such a case, the term corresponding to the 
destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_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 the four neighbouring 
$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 (\S\ref{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 
(see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical 
diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} 
namelist parameter to true and defining the \key{zdftke} CPP key all together.

The KPP 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{np\_zdfevd} should not be used with the KPP 
scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module...

% ================================================================
% Double Diffusion Mixing
% ================================================================
\section  [Double Diffusion Mixing (\textit{zdfddm} - \key{zdfddm})]
      {Double Diffusion Mixing (\mdl{zdfddm} module - \key{zdfddm})}
\label{ZDF_ddm}

%-------------------------------------------namddm--------------------------------------------------------
\namdisplay{namddm}
%--------------------------------------------------------------------------------------------------------------

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{Merryfield1999} 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_zdfddm_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 \S\ref{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_zdfddm_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_zdfddm_f_T}
A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
\end{align}

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \label{Fig_zdfddm}  \begin{center}
\includegraphics[width=0.99\textwidth]{./TexFiles/Figures/Fig_zdfddm.pdf}
\caption {From \citet{Merryfield1999} : (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. }
\end{center}    \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

The factor 0.7 in \eqref{Eq_zdfddm_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 ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following  \citet{Merryfield1999}, 
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_zdfddm_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}   
\\          \label{Eq_zdfddm_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 \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ 
are illustrated in Fig.~\ref{Fig_zdfddm}. 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).

% ================================================================
% Bottom Friction
% ================================================================
\section  [Bottom Friction (\textit{zdfbfr})]
      {Bottom Friction (\mdl{zdfbfr} module)}
\label{ZDF_bfr}

%--------------------------------------------nambfr--------------------------------------------------------
\namdisplay{nambfr}
%--------------------------------------------------------------------------------------------------------------

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:
\begin{equation} \label{Eq_zdfbfr_flux}
A^{vm} \left( \partial \textbf{U}_h / \partial z \right) = \textbf{F}_h
\end{equation}
where $\textbf{F}_h$ is supposed to represent the horizontal momentum flux 
outside the logarithmic turbulent boundary layer (thickness of the order of 
1~m in the ocean). How $\textbf{F}_h$ 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_zdfbfr_flux2}
\frac{\partial u \; (k)}{\partial t} = \frac{1}{e_{3u}} \left[ A^{vm} \; (k) \frac{U \; (k-1) - U \; (k)}{e_{3uw} \; (k-1)} - F_u \right] \approx - \frac{F_u}{e_{3u}}
\end{equation}
For example, 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 $\textbf{F}_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.

% -------------------------------------------------------------------------------------------------------------
%       Linear Bottom Friction
% -------------------------------------------------------------------------------------------------------------
\subsection{Linear Bottom Friction (\np{nbotfr} = 1) }
\label{ZDF_bfr_linear}

The linear bottom friction parameterisation assumes that the bottom friction 
is proportional to the interior velocity (i.e. the velocity of the last model level):
\begin{equation} \label{Eq_zdfbfr_linear}
\textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
\end{equation}
where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean 
layer and $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{Weatherly1984}. 
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{Gill1982}, 
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{bfric1} (namelist parameter).

In the code, the bottom friction is imposed by updating the value of the 
vertical eddy coefficient at the bottom level. Indeed, the discrete formulation 
of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that 
$\textbf {U}_h =0$ below the ocean floor, leads to
\begin{equation} \label{Eq_zdfbfr_linKz}
\begin{split}
A_u^{vm} &= r\;e_{3uw}\\
A_v^{vm} &= r\;e_{3vw}\\
\end{split}
\end{equation}

This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$ 
used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and 
leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets 
$r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background 
vertical eddy coefficient, and a no-slip boundary condition is imposed. 
Note that this latter choice generally leads to an underestimation of the 
bottom friction: for example with a deepest level thickness of $200~m$ 
and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient 
is only $r=10^{-6}$m.s$^{-1}$.

% -------------------------------------------------------------------------------------------------------------
%       Non-Linear Bottom Friction
% -------------------------------------------------------------------------------------------------------------
\subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)}
\label{ZDF_bfr_nonlinear}

The non-linear bottom friction parameterisation assumes that the bottom 
friction is quadratic: 
\begin{equation} \label{Eq_zdfbfr_nonlinear}
\textbf {F}_h = \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 
\end{equation}

with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$ 
the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient, 
and $e_b $ a 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{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$, 
while the FRAM experiment \citep{Killworth1992} uses $e_b =0$ 
and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been 
set as default values (\np{bfric2} and \np{bfeb2} namelist parameters).

As for the linear case, the bottom friction is imposed in the code by 
updating the value of the vertical eddy coefficient at the bottom level:
\begin{equation} \label{Eq_zdfbfr_nonlinKz}
\begin{split}
A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^
{1/2}\\
A_v^{vm} &=C_D\; e_{3uw} \left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
\end{split}
\end{equation}

This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the 
non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, 
and $e_b$ =\np{bfeb2}.

% ================================================================