source: trunk/NEMO/DOC/BETA/Chapters/Chap_ZDF.tex @ 705

% ================================================================
% 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.

% ================================================================
% 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 parameterization 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 \mdl{zdfini} module and performed in \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 \mdl{dynzdf} and \mdl{trazdf} modules, respectively. These trends can be computed using either a forward time scheme (cpp variable 
\np{np\_zdfexp} or a backward time scheme (default option) depending on the magnitude of the mixing coefficients used, and thus of the formulation used (see \S\ref{DOM_nxt}).

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

When the \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 to use this option only in process 
studies, not in basin scale simulation. 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 \np{avm0} and \np{avt0} namelist parameters. In any case, do not use 
values smaller that those associated to 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 of the vertical momentum and tracer eddy coefficients is set. The vertical mixing coefficients are diagnosed from the large scale variables computed by the model (order 0.5 closure scheme). \textit{In situ} measurements allow 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 turbulent eddy coefficients and the local Richardson number ($i.e.$ ratio of stratification over 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 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 coefficients can be modified by setting \np{avmri}, \np{alp} and \np{nric} namelist parameter, respectively.

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

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

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 atmospheric cases, adapted by \citet{Gaspar1990} for oceanic cases, 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$ designates the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}), 
$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing turbulent 
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 parameter \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 be optionally applied to $R_i$. Nevertheless it is an obsolescent option that is notrecommanded.  The choice of $P_{rt} $ is controlled by \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 \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 obtained in a stable stratified region with constant values of the brunt-Vais\"{a}l\"{a} frequency. The resulting turbulent 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 has no sense for local unstable 
stratification and the computation no longer uses the whole information contained in the vertical density profile. To overcome this drawbacks, \citet{Madec1998} introduces the \np{nmxl}=2 or 3 cases, which add of an hypothesis on 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 imposed \eqref{Eq_tke_mxl_constraint} constraint, we introduce two additonnal length scal: $l_{up}$ and $l_{dwn}$, the upward and downward length scale, and evaluate the dissipation and mixing turbulent length scales as (caution here we use the numerical indexation):
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \label{Fig_mixing_length}  \begin{center}
\includegraphics[width=1.00\textwidth]{./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 compute 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 turbulent length scales take a same value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the \np{nmxl}=2 case, the dissipation and mixing turbulent 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 bottom value 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 the 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 
\citet{Gargett1984} one 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 KKP scheme has been implemented by J. Chanut ...
\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: 
they must be parameterized. Three parameterisations are available to deal 
with convective processes: either 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}=T) }
\label{ZDF_npc}

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


%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!ht] \label{Fig_npc}    \begin{center}
\includegraphics[width=0.90\textwidth]{./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 algorithm is used when \np{ln\_zdfnpc}=T. 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}. This algorithm is an iterative process used in the following way 
(Fig. \ref{Fig_npc}): going from the top of the ocean towards the bottom, the first 
instability is searched. Assume in the following that the instability is 
located between levels $k$ and $k+1$. The two levels are vertically mixed, for 
potential temperature and salinity, 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 two by two statically unstable levels. The latter procedure cannot converge with a finite number 
of iterations for some vertical profiles while the algorithm used in OPA 
converges for any profile in a number of iterations 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 studying the deep water formation in 
the north-western Mediterranean Sea \citep{Madec1991a, Madec1991b, Madec1991c}.

Note that in this algorithm the potential density referenced to the sea 
surface is used to check whether the density profile is stable or not. 
Moreover, the mixing in 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. Moreover, temperature and 
salinity, and thus density, are mixed, but the corresponding velocity fields 
remain unchanged. When using a Richardson 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 algorithm is used with constant 
vertical eddy viscosity, spurious solution can occur as the vertical 
momentum diffusion remains small even after a static adjustment. In that 
latter case, we recommend to add momentum mixing in a manner that mimics the 
mixing in temperature and salinity \citep{Speich1992, Speich1996}.

%AMT This presentation fails to mention the big drawback of this scheme: many water masses of the world ocean, especially AABW, are unstable when represented in surface-referenced potential density sigma_0. The scheme erroneously mixes them up.

% -------------------------------------------------------------------------------------------------------------
%       Enhanced Vertical Diffusion 
% -------------------------------------------------------------------------------------------------------------
\subsection{Enhanced Vertical Diffusion (\np{ln\_zdfevd}=T)}
\label{ZDF_evd}

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

The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd} is 
defined. In this case, the vertical eddy mixing coefficients are assigned to be very large (a typical value is $1\;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) mixing coefficients.

In practice, when $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$ are set to a large value, 
\np{avevd}, and  if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm} $ . Typical value for $avevd$ is in-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 case of static instabilities. It requires the use of an implicit time stepping on vertical diffusion terms (i.e. \np{ln\_zdfexp}=F). 

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

The turbulent closure scheme presented in \S\ref{ZDF_tke} and used when the 
\key{zdftke} is defined allows, in theory, to deal with 
statically unstable density profiles. In such a 
case, the term of destruction of turbulent kinetic energy through 
stratification in \eqref{Eq_zdftke_e} becomes a source term as $N^2$ is negative. It 
results large values of both $A_T^{vT}$ and the four neighbouring$A_u^{vm} 
{and}\;A_v^{vm} $ (up to $1\;m^2s^{-1})$ that are able to restore the 
static stability of the water column in a way similar to that of the 
enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). Nevertheless, the 
eddy coefficients computed by the turbulent scheme do usually not exceed
$10^{-2}m.s^{-1}$ in the vicinity of the sea surface (first ocean layer) due to 
the bound of the turbulent length scale by the distance to the sea surface 
(see {\S}VI.7-c). It can thus be useful to combine the enhanced vertical 
diffusion with the turbulent closure, i.e. defining \np{np\_zdfevd} and 
\key{zdftke} CPP variables all together.

The KPP scheme 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


% ================================================================
% Double Diffusion Mixing
% ================================================================
\section{Double Diffusion Mixing (\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 oceans.  \citet{Merryfield1999} include a parameterization 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 buoyancy ratio $R_\rho = \alpha \partial_z T / \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}   
\\
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]{./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^{vT}$ for temperature and salt in regions of diffusive convection. Heavy curves denote the Federov parameterization and thin curves the Kelley parameterization. The latter is not implemented in \NEMO }
\end{center}    \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

The factor 0.7 in \eqref{Eq_zdfddm_f} reflects the measured ratio $\alpha F_T /\beta F_S \approx  0.7$ of buoyancy fluxes due to transport of heat and salt (e.g., McDougall and Taylor 1984). 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}   
\\
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 dependences of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d} on $R_\rho$ are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing $R_\rho$ at each 
grid point and 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}
\label{ZDF_bfr}

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

Both 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 OPA: 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 OPA.

% -------------------------------------------------------------------------------------------------------------
%       Linear Bottom Friction
% -------------------------------------------------------------------------------------------------------------
\subsection{Linear Bottom Friction}
\label{ZDF_bfr_linear}

%-------------------------------------------nambfr--------------------------------------------------------
\begin{flushright}
(\textbf{namelist} !nbotfr :\textit{ nbotfr = 0, = 1 or = 3})
\end{flushright}

The linear bottom friction parameterization 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$ a friction coefficient expressed in m.s$^{-1}$. This 
coefficient is generally estimated by setting a typical decay time $\tau $ in the 
deep ocean, $r = H / \tau$. Commonly accepted values of $\tau$ are of the 
order of 100 to 200 days \citep{Weatherly1984}. A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ 
corresponding 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). With a drag coefficient $C_D = 0.002$, a typical value of tidal currents $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 
OPA. 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 specified 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$ inside the bottom, leads to
\begin{equation} \label{Eq_zdfbfr_linKz}
\begin{split}
A_u^{vm} &= r\;e_{3uw}\\
A_v^{vm} &= r\;e_{3uw}\\
\end{split}
\end{equation}

Such an update is done in \mdl{zdfbfr} when \np{nbotfr}=1 and the value of $r$ used is 
\np{bfric1}. Setting \np{nbotfr}=3 is equivalent to set $r=0$ and leads to a free-slip bottom boundary condition, 
while setting \np{nbotfr}=0 imposes $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the 
background vertical eddy coefficient: a no-slip boundary condition is used. 
Note that this latter choice generally leads to an underestimation of the 
bottom friction: for 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}
\label{ZDF_bfr_nonlinear}

\begin{center}
(\textbf{namelist} !nbotfr : \textit{nbotfr = 2})
\end{center}

The non-linear bottom friction parameterization 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 value (\np{bfric2} and \np{bfeb2} namelist parameters).

As for the linear case, the bottom friction is specified 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+1,j}\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}).

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