New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 994 for trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex – NEMO

Ignore:
Timestamp:
2008-05-28T11:01:09+02:00 (16 years ago)
Author:
gm
Message:

trunk - add steven correction + several other things + rename BETA into TexFiles?

Location:
trunk/DOC/TexFiles
Files:
1 copied
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex

    r817 r994  
    77 
    88%gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. 
     9\gmcomment{Steven remark : 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! 
     10Gurvan :  I kept "turbulent closure"...} 
     11\gmcomment{Steven bis : parameterization is the american spelling, parameterisation is the british} 
     12 
    913 
    1014% ================================================================ 
    … …  
    1418\label{ZDF_zdf} 
    1519 
    16 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}). 
     20The discrete form of the ocean subgrid scale physics has been presented in  
     21\S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries,  
     22the turbulent fluxes of momentum, heat and salt have to be defined. At the  
     23surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}),  
     24while at the bottom they are set to zero for heat and salt, unless a geothermal  
     25flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl}  
     26defined, see \S\ref{TRA_bbc}), and specified through a bottom friction  
     27parameterization for momentum (see \S\ref{ZDF_bfr}). 
    1728 
    1829In this section we briefly discuss the various choices offered to compute  
    19 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,  
    20 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  
    21 \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}). 
     30the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ ,  
     31$A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$-  
     32points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These  
     33coefficients can be assumed to be either constant, or a function of the local  
     34Richardson number, or computed from a turbulent closure model (either  
     35TKE or KPP formulation). The computation of these coefficients is initialized  
     36in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or  
     37\mdl{zdfkpp} modules. The trends due to the vertical momentum and tracer  
     38diffusion, including the surface forcing, are computed and added to the  
     39general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.  
     40These trends can be computed using either a forward time stepping scheme  
     41(namelist parameter \np{np\_zdfexp}=true) or a backward time stepping  
     42scheme (\np{np\_zdfexp}=false) depending on the magnitude of the mixing  
     43coefficients, and thus of the formulation used (see \S\ref{DOM_nxt}). 
    2244 
    2345% ------------------------------------------------------------------------------------------------------------- 
    … …  
    3052%-------------------------------------------------------------------------------------------------------------- 
    3153 
    32 When the \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to  
    33 constant values over the whole ocean. This is the crudest way to define the  
    34 vertical ocean physics. It is recommended to use this option only in process  
    35 studies, not in basin scale simulation. Typical values used in this case  
    36 are: 
     54When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients  
     55are set to constant values over the whole ocean. This is the crudest way to define  
     56the vertical ocean physics. It is recommended that this option is only used in  
     57process studies, not in basin scale simulations. Typical values used in this case are: 
    3758\begin{align*}  
    3859A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}  \\ 
    … …  
    4162\end{align*} 
    4263 
    43 These values are set through \np{avm0} and \np{avt0} namelist parameters. In any case, do not use  
    44 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. 
     64These values are set through the \np{avm0} and \np{avt0} namelist parameters.  
     65In all cases, do not use values smaller that those associated with the molecular  
     66viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum,  
     67$\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity. 
    4568 
    4669 
    … …  
    5578%-------------------------------------------------------------------------------------------------------------- 
    5679 
    57 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: 
     80When \key{zdfric} is defined, a local Richardson number dependent formulation  
     81for the vertical momentum and tracer eddy coefficients is set. The vertical mixing  
     82coefficients are diagnosed from the large scale variables computed by the model.  
     83\textit{In situ} measurements have been used to link vertical turbulent activity to  
     84large scale ocean structures. The hypothesis of a mixing mainly maintained by the  
     85growth of Kelvin-Helmholtz like instabilities leads to a dependency between the  
     86vertical turbulence eddy coefficients and the local Richardson number ($i.e.$ the  
     87ratio of stratification to vertical shear). Following \citet{PacPhil1981}, the following  
     88formulation has been implemented: 
    5889\begin{equation} \label{Eq_zdfric} 
    5990   \left\{      \begin{aligned} 
    … …  
    6394   \end{aligned}    \right. 
    6495\end{equation} 
    65 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. 
     96where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson  
     97number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),  
     98$A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the  
     99constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$  
     100is the maximum value that can be reached by the coefficient when $Ri\leq 0$,  
     101$a=5$ and $n=2$. The last three values can be modified by setting the  
     102\np{avmri}, \np{alp} and \np{nric} namelist parameters, respectively. 
    66103 
    67104% ------------------------------------------------------------------------------------------------------------- 
    … …  
    75112%-------------------------------------------------------------------------------------------------------------- 
    76113 
    77 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: 
     114The vertical eddy viscosity and diffusivity coefficients are computed from a TKE  
     115turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent  
     116kinetic energy, and a closure assumption for the turbulence length scales. This  
     117turbulent closure model has been developed by \citet{Bougeault1989} in the  
     118atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and  
     119embedded in OPA by \citet{Blanke1993} for equatorial Atlantic simulations. Since  
     120then, significant modifications have been introduced by \citet{Madec1998} in both  
     121the implementation and the formulation of the mixing length scale. The time  
     122evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical  
     123shear, its destruction through stratification, its vertical diffusion, and its dissipation  
     124of \citet{Kolmogorov1942} type: 
    78125\begin{equation} \label{Eq_zdftke_e} 
    79126\frac{\partial \bar{e}}{\partial t} =  
    … …  
    87134\begin{equation} \label{Eq_zdftke_kz} 
    88135   \begin{split} 
    89          A^{vm} &= C_k\  l_k\  \sqrt {\bar{e}} 
    90     \\ 
     136         A^{vm} &= C_k\  l_k\  \sqrt {\bar{e}}     \\ 
    91137         A^{vT} &= A^{vm} / P_{rt} 
    92138   \end{split} 
    93139\end{equation} 
    94 where $N$ designates the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),  
    95 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing turbulent  
    96 length scales, $P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and  
    97 $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  
    98 set to unity or, following \citet{Blanke1993}, be a function of the  
    99 local Richardson number, $R_i $: 
     140where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),  
     141$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,  
     142$P_{rt} $ is the Prandtl number. The constants $C_k = \sqrt {2} /2$ and  
     143$C_\epsilon = 0.1$ are designed to deal with vertical mixing at any depth  
     144\citep{Gaspar1990}. They are set through namelist parameters \np{ediff}  
     145and \np{ediss}. $P_{rt} $ can be set to unity or, following \citet{Blanke1993},  
     146be a function of the local Richardson number, $R_i $: 
    100147\begin{align*} \label{Eq_prt} 
    101148P_{rt} = \begin{cases} 
    … …  
    105152            \end{cases} 
    106153\end{align*} 
    107 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. 
    108  
    109 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}: 
     154Note that a horizontal Shapiro filter can optionally be applied to $R_i$.  
     155However it is an obsolescent option that is not recommended.   
     156The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter. 
     157 
     158For computational efficiency, the original formulation of the turbulence length  
     159scales proposed by \citet{Gaspar1990} has been simplified. Four formulations  
     160are proposed, the choice of which is controlled by the \np{nmxl} namelist  
     161parameter. The first two are based on the following first order approximation  
     162\citep{Blanke1993}: 
    110163\begin{equation} \label{Eq_tke_mxl0_1} 
    111164l_k = l_\epsilon = \sqrt {2 \bar e} / N 
    112165\end{equation} 
    113 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  
    114 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: 
     166which is valid in a stable stratified region with constant values of the brunt- 
     167Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance  
     168to 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  
     169drawbacks: it makes no sense for locally unstable stratification and the  
     170computation no longer uses all the information contained in the vertical density  
     171profile. To overcome these drawbacks, \citet{Madec1998} introduces the  
     172\np{nmxl}=2 or 3 cases, which add an extra assumption concerning the vertical  
     173gradient of the computed length scale. So, the length scales are first evaluated  
     174as in \eqref{Eq_tke_mxl0_1} and then bounded such that: 
    115175\begin{equation} \label{Eq_tke_mxl_constraint} 
    116176\frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1 
    117177\qquad \text{with }\  l =  l_k = l_\epsilon 
    118178\end{equation} 
    119  
    120 \eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length scale cannot be  
    121 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  
    122 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): 
     179\eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length  
     180scale cannot be larger than the variations of depth. It provides a better  
     181approximation of the \citet{Gaspar1990} formulation while being much less  
     182time consuming. In particular, it allows the length scale to be limited not only  
     183by the distance to the surface or to the ocean bottom but also by the distance  
     184to a strongly stratified portion of the water column such as the thermocline  
     185(Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint}  
     186constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$,  
     187the upward and downward length scales, and evaluate the dissipation and  
     188mixing turbulence length scales as (and note that here we use numerical  
     189indexing): 
    123190%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    124191\begin{figure}[!t] \label{Fig_mixing_length}  \begin{center} 
    … …  
    130197\begin{equation} \label{Eq_tke_mxl2} 
    131198\begin{aligned} 
    132  l_{up  }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)} \right) 
     199 l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3T}^{(k)}\ \ \ \; \right) 
    133200    \quad &\text{ from $k=1$ to $jpk$ }\ \\ 
    134201 l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3T}^{(k-1)}  \right)    
    … …  
    136203\end{aligned} 
    137204\end{equation} 
    138 where $l^{(k)}$ is compute using \eqref{Eq_tke_mxl0_1}, $i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$. 
    139  
    140 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}: 
     205where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1},  
     206$i.e.$ $l^{(k)} = \sqrt {2 \bar e^{(k)} / N^{(k)} }$. 
     207 
     208In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same  
     209value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the  
     210\np{nmxl}=2 case, the dissipation and mixing turbulence length scales are give  
     211as in \citet{Gaspar1990}: 
    141212\begin{equation} \label{Eq_tke_mxl_gaspar} 
    142213\begin{aligned} 
    … …  
    149220stress field: $\bar{e}=ebb\;\left| \tau \right|$ ($ebb=60$ by default)  
    150221with a minimal threshold of $emin0=10^{-4}~m^2.s^{-2}$ (namelist 
    151 parameters). Its bottom value is assumed to be equal to the value of the level just above. The time  
    152 integration of the $\bar{e}$ equation may formally lead to negative values  
    153 because the numerical scheme does not ensure the positivity. To overcome  
    154 this problem, a cut-off in the minimum value of $\bar{e}$ is used. Following  
    155 \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  
    156 \citet{Gargett1984} one for the diffusion in the thermocline and deep ocean  
    157 $(A^{vT} = 10^{-3} / N)$. In addition, a  
    158 cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical  
     222parameters). Its value at the bottom of the ocean is assumed to be  
     223equal to the value of the level just above. The time integration of the  
     224$\bar{e}$ equation may formally lead to negative values because the  
     225numerical scheme does not ensure its positivity. To overcome this  
     226problem, a cut-off in the minimum value of $\bar{e}$ is used. Following  
     227\citet{Gaspar1990}, the cut-off value is set to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$.  
     228This allows the subsequent formulations to match that of\citet{Gargett1984}  
     229for the diffusion in the thermocline and deep ocean :  $(A^{vT} = 10^{-3} / N)$.  
     230In addition, a cut-off is applied on $A^{vm}$ and $A^{vT}$ to avoid numerical  
    159231instabilities associated with too weak vertical diffusion. They must be  
    160232specified at least larger than the molecular values, and are set through  
    … …  
    172244 
    173245The KKP scheme has been implemented by J. Chanut ... 
     246 
    174247\colorbox{yellow}{Add a description of KPP here.} 
    175248 
    … …  
    188261occur at particular ocean grid points. In nature, convective processes  
    189262quickly re-establish the static stability of the water column. These  
    190 processes have been removed from the model via the hydrostatic assumption:  
    191 they must be parameterized. Three parameterisations are available to deal  
    192 with convective processes: either a non-penetrative convective adjustment or  
    193 an enhanced vertical diffusion, or/and the use of a turbulent closure  
    194 scheme. 
     263processes have been removed from the model via the hydrostatic  
     264assumption so they must be parameterized. Three parameterizations  
     265are available to deal with convective processes: a non-penetrative  
     266convective adjustment or an enhanced vertical diffusion, or/and the  
     267use of a turbulent closure scheme. 
    195268 
    196269% ------------------------------------------------------------------------------------------------------------- 
    … …  
    207280 
    208281%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    209 \begin{figure}[!ht] \label{Fig_npc}    \begin{center} 
     282\begin{figure}[!htb] \label{Fig_npc}   \begin{center} 
    210283\includegraphics[width=0.90\textwidth]{./Figures/Fig_npc.pdf} 
    211 \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.} 
     284\caption {Example of an unstable density profile treated by the non penetrative  
     285convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from  
     286the surface to the bottom. It is found to be unstable between levels 3 and 4.  
     287They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5  
     288are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are  
     289mixed. The $1^{st}$ step ends since the density profile is then stable below  
     290the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same  
     291procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile  
     292is checked. It is found stable: end of algorithm.} 
    212293\end{center}   \end{figure} 
    213294%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    214295 
    215 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  
    216 (Fig. \ref{Fig_npc}): going from the top of the ocean towards the bottom, the first  
    217 instability is searched. Assume in the following that the instability is  
    218 located between levels $k$ and $k+1$. The two levels are vertically mixed, for  
    219 potential temperature and salinity, conserving the heat and salt contents of  
    220 the water column. The new density is then computed by a linear  
    221 approximation. If the new density profile is still unstable between levels  
    222 $k+1$ and $k+2$, levels $k$, $k+1$ and $k+2$ are then mixed. This process is repeated until  
    223 stability is established below the level $k$ (the mixing process can go down to  
    224 the ocean bottom). The algorithm is repeated to check if the density profile  
     296The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}=true.  
     297It is applied at each \np{nnpc1} time step and mixes downwards instantaneously  
     298the statically unstable portion of the water column, but only until the density  
     299structure becomes neutrally stable ($i.e.$ until the mixed portion of the water  
     300column has \textit{exactly} the density of the water just below) \citep{Madec1991a}.  
     301The associated algorithm is an iterative process used in the following way  
     302(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is  
     303found. Assume in the following that the instability is located between levels  
     304$k$ and $k+1$. The potential temperature and salinity in the two levels are  
     305vertically mixed, conserving the heat and salt contents of the water column.  
     306The new density is then computed by a linear approximation. If the new  
     307density profile is still unstable between levels $k+1$ and $k+2$, levels $k$,  
     308$k+1$ and $k+2$ are then mixed. This process is repeated until stability is  
     309established below the level $k$ (the mixing process can go down to the  
     310ocean bottom). The algorithm is repeated to check if the density profile  
    225311between level $k-1$ and $k$ is unstable and/or if there is no deeper instability. 
    226312 
    227 This algorithm is significantly different from mixing two by two statically unstable levels. The latter procedure cannot converge with a finite number  
    228 of iterations for some vertical profiles while the algorithm used in OPA  
    229 converges for any profile in a number of iterations less than the number of  
    230 vertical levels. This property is of paramount importance as pointed out by  
    231 \citet{Killworth1989}: it avoids the existence of permanent and unrealistic  
    232 static instabilities at the sea surface. This non-penetrative convective  
    233 algorithm has been proved successful in studying the deep water formation in  
    234 the north-western Mediterranean Sea \citep{Madec1991a, Madec1991b, Madec1991c}. 
    235  
    236 Note that in this algorithm the potential density referenced to the sea  
    237 surface is used to check whether the density profile is stable or not.  
    238 Moreover, the mixing in potential density is assumed to be linear. This  
    239 assures the convergence of the algorithm even when the equation of state is  
    240 non-linear. Small static instabilities can thus persist due to cabbeling:  
    241 they will be treated at the next time step. Moreover, temperature and  
    242 salinity, and thus density, are mixed, but the corresponding velocity fields  
    243 remain unchanged. When using a Richardson dependent eddy viscosity, the  
    244 mixing of momentum is done through the vertical diffusion: after a static  
    245 adjustment, the Richardson number is zero and thus the eddy viscosity  
    246 coefficient is at a maximum. When this algorithm is used with constant  
    247 vertical eddy viscosity, spurious solution can occur as the vertical  
    248 momentum diffusion remains small even after a static adjustment. In that  
    249 latter case, we recommend to add momentum mixing in a manner that mimics the  
    250 mixing in temperature and salinity \citep{Speich1992, Speich1996}. 
    251  
    252 %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. 
     313This algorithm is significantly different from mixing statically unstable levels  
     314two by two. The latter procedure cannot converge with a finite number  
     315of iterations for some vertical profiles while the algorithm used in \NEMO  
     316converges for any profile in a number of iterations which is less than the  
     317number of vertical levels. This property is of paramount importance as  
     318pointed out by \citet{Killworth1989}: it avoids the existence of permanent  
     319and unrealistic static instabilities at the sea surface. This non-penetrative  
     320convective algorithm has been proved successful in studies of the deep  
     321water formation in the north-western Mediterranean Sea  
     322\citep{Madec1991a, Madec1991b, Madec1991c}. 
     323 
     324Note that in the current implementation of this algorithm presents several  
     325limitations. First, potential density referenced to the sea surface is used to  
     326check whether the density profile is stable or not. This is a strong  
     327simplification which leads to large errors for realistic ocean simulations.  
     328Indeed, many water masses of the world ocean, especially Antarctic Bottom 
     329Water, are unstable when represented in surface-referenced potential density.  
     330The scheme will erroneously mix them up. Second, the mixing of potential  
     331density is assumed to be linear. This assures the convergence of the algorithm  
     332even when the equation of state is non-linear. Small static instabilities can thus  
     333persist due to cabbeling: they will be treated at the next time step.  
     334Third, temperature and salinity, and thus density, are mixed, but the  
     335corresponding velocity fields remain unchanged. When using a Richardson  
     336Number dependent eddy viscosity, the mixing of momentum is done through  
     337the vertical diffusion: after a static adjustment, the Richardson Number is zero  
     338and thus the eddy viscosity coefficient is at a maximum. When this convective  
     339adjustment algorithm is used with constant vertical eddy viscosity, spurious  
     340solutions can occur since the vertical momentum diffusion remains small even  
     341after a static adjustment. In that case, we recommend the addition of momentum  
     342mixing in a manner that mimics the mixing in temperature and salinity  
     343\citep{Speich1992, Speich1996}. 
    253344 
    254345% ------------------------------------------------------------------------------------------------------------- 
    … …  
    263354%-------------------------------------------------------------------------------------------------------------- 
    264355 
    265 The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd} is  
    266 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  
    267 tracers (\np{n\_evdm}=1) mixing coefficients. 
    268  
    269 In practice, when $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$ are set to a large value,  
    270 \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).  
     356The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd}=true.  
     357In this case, the vertical eddy mixing coefficients are assigned very large values  
     358(a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable  
     359($i.e.$ when the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{Lazar1997, Lazar1999}.  
     360This is done either on tracers only (\np{n\_evdm}=0) or on both momentum and  
     361tracers (\np{n\_evdm}=1). 
     362 
     363In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and  
     364if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$  
     365values also, are set equal to the namelist parameter \np{avevd}. A typical value  
     366for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterization of  
     367convective processes is less time consuming than the convective adjustment  
     368algorithm presented above when mixing both tracers and momentum in the  
     369case of static instabilities. It requires the use of an implicit time stepping on  
     370vertical diffusion terms (i.e. \np{ln\_zdfexp}=false).  
    271371 
    272372% ------------------------------------------------------------------------------------------------------------- 
    … …  
    276376\label{ZDF_tcs} 
    277377 
    278 The turbulent closure scheme presented in \S\ref{ZDF_tke} and used when the  
    279 \key{zdftke} is defined allows, in theory, to deal with  
    280 statically unstable density profiles. In such a  
    281 case, the term of destruction of turbulent kinetic energy through  
    282 stratification in \eqref{Eq_zdftke_e} becomes a source term as $N^2$ is negative. It  
    283 results large values of both $A_T^{vT}$ and the four neighbouring$A_u^{vm}  
    284 {and}\;A_v^{vm} $ (up to $1\;m^2s^{-1})$ that are able to restore the  
    285 static stability of the water column in a way similar to that of the  
    286 enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). Nevertheless, the  
    287 eddy coefficients computed by the turbulent scheme do usually not exceed 
    288 $10^{-2}m.s^{-1}$ in the vicinity of the sea surface (first ocean layer) due to  
    289 the bound of the turbulent length scale by the distance to the sea surface  
    290 (see {\S}VI.7-c). It can thus be useful to combine the enhanced vertical  
    291 diffusion with the turbulent closure, i.e. defining \np{np\_zdfevd} and  
    292 \key{zdftke} CPP variables all together. 
    293  
    294 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 
    295  
     378The TKE turbulent closure scheme presented in \S\ref{ZDF_tke} and used  
     379when the \key{zdftke} is defined, in theory solves the problem of statically  
     380unstable density profiles. In such a case, the term corresponding to the  
     381destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e}  
     382becomes a source term, since $N^2$ is negative. It results in large values of  
     383$A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring  
     384$A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values  
     385restore the static stability of the water column in a way similar to that of the  
     386enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). However,  
     387in the vicinity of the sea surface (first ocean layer), the eddy coefficients  
     388computed by the turbulence scheme do not usually exceed $10^{-2}m.s^{-1}$,  
     389because the mixing length scale is bounded by the distance to the sea surface  
     390(see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical  
     391diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc}  
     392namelist parameter to true and defining the \key{zdftke} CPP key all together. 
     393 
     394The KPP turbulent closure scheme already includes enhanced vertical diffusion  
     395in the case of convection, as governed by the variables $bvsqcon$ and $difcon$  
     396found in \mdl{zdfkpp}, therefore \np{np\_zdfevd} should not be used with the KPP  
     397scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module... 
    296398 
    297399% ================================================================ 
    … …  
    306408%-------------------------------------------------------------------------------------------------------------- 
    307409 
    308 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.  
     410Double diffusion occurs when relatively warm, salty water overlies cooler, fresher  
     411water, or vice versa. The former condition leads to salt fingering and the latter  
     412to diffusive convection. Double-diffusive phenomena contribute to diapycnal  
     413mixing in extensive regions of the ocean.  \citet{Merryfield1999} include a  
     414parameterization of such phenomena in a global ocean model and show that  
     415it leads to relatively minor changes in circulation but exerts significant regional  
     416influences on temperature and salinity.  
    309417 
    310418Diapycnal mixing of S and T are described by diapycnal diffusion coefficients  
    … …  
    313421    &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS} 
    314422\end{align*} 
    315 where subscript $f$ represents mixing by salt fingering,  
    316 $d$ by diffusive convection, and $o$ by processes other than  
    317 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): 
     423where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,  
     424and $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$,  
     425where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline  
     426contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt  
     427fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981): 
    318428\begin{align} \label{Eq_zdfddm_f} 
    319429A_f^{vS} &=    \begin{cases} 
    … …  
    321431   0                              &\text{otherwise}  
    322432            \end{cases}    
    323 \\ 
     433\\           \label{Eq_zdfddm_f_T} 
    324434A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho  
    325435\end{align} 
    … …  
    328438\begin{figure}[!t] \label{Fig_zdfddm}  \begin{center} 
    329439\includegraphics[width=0.99\textwidth]{./Figures/Fig_zdfddm.pdf} 
    330 \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 } 
     440\caption {From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$  
     441and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy  
     442curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves  
     443$A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and  
     444$A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy  
     445curves denote the Federov parameterization and thin curves the Kelley  
     446parameterization. The latter is not implemented in \NEMO. } 
    331447\end{center}    \end{figure} 
    332448%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    333449 
    334 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}$. 
     450The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio  
     451$\alpha F_T /\beta F_S \approx  0.7$ of buoyancy flux of heat to buoyancy  
     452flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following  \citet{Merryfield1999},  
     453we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$. 
    335454 
    336455To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used:  
    337 \begin{align} \label{Eq_zdfddm_d} 
     456\begin{align}  \label{Eq_zdfddm_d} 
    338457A_d^{vT} &=    \begin{cases} 
    339458   1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)} 
    … …  
    341460   0                       &\text{otherwise}  
    342461            \end{cases}    
    343 \\ 
     462\\          \label{Eq_zdfddm_d_S} 
    344463A_d^{vS} &=    \begin{cases} 
    345464   A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) 
    … …  
    351470\end{align} 
    352471 
    353 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  
    354 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). 
    355  
     472The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$  
     473are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing  
     474$R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the  
     475same time as $N^2$ is computed. This avoids duplication in the computation of  
     476$\alpha$ and $\beta$ (which is usually quite expensive). 
    356477 
    357478% ================================================================ 
    … …  
    366487%-------------------------------------------------------------------------------------------------------------- 
    367488 
    368 Both surface momentum flux (wind stress) and the bottom momentum flux  
    369 (bottom friction) enter the equations as a condition on the vertical  
     489Both the surface momentum flux (wind stress) and the bottom momentum  
     490flux (bottom friction) enter the equations as a condition on the vertical  
    370491diffusive flux. For the bottom boundary layer, one has: 
    371492\begin{equation} \label{Eq_zdfbfr_flux} 
    … …  
    3764971~m in the ocean). How $\textbf{F}_h$ influences the interior depends on the  
    377498vertical resolution of the model near the bottom relative to the Ekman layer  
    378 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  
    379 frequency $f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient  
    380 $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  
    381 viscous forces (either wind stress or bottom friction) as a body force over  
    382 the depth of the top or bottom model layer. To illustrate this, consider the  
    383 equation for $u$ at $k$, the last ocean level: 
     499depth. For example, in order to obtain an Ekman layer depth  
     500$d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient  
     501$A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency  
     502$f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient  
     503$A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.  
     504When the vertical mixing coefficient is this small, using a flux condition is  
     505equivalent to entering the viscous forces (either wind stress or bottom friction)  
     506as a body force over the depth of the top or bottom model layer. To illustrate  
     507this, consider the equation for $u$ at $k$, the last ocean level: 
    384508\begin{equation} \label{Eq_zdfbfr_flux2} 
    385509\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}} 
    386510\end{equation} 
    387 For example, if the bottom layer thickness is 200~m, the Ekman transport will be  
    388 distributed over that depth. On the other hand, if the vertical resolution  
     511For example, if the bottom layer thickness is 200~m, the Ekman transport will  
     512be distributed over that depth. On the other hand, if the vertical resolution  
    389513is high (1~m or less) and a turbulent closure model is used, the turbulent  
    390514Ekman layer will be represented explicitly by the model. However, the  
    391515logarithmic layer is never represented in current primitive equation model  
    392 applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $. Two  
    393 choices are available in OPA: a linear and a quadratic bottom friction. Note  
    394 that in both cases, the rotation between the interior velocity and the  
    395 bottom friction is neglected in the present release of OPA. 
     516applications: it is \emph{necessary} to parameterize the flux $\textbf{F}_h $.  
     517Two choices are available in \NEMO: a linear and a quadratic bottom friction.  
     518Note that in both cases, the rotation between the interior velocity and the  
     519bottom friction is neglected in the present release of \NEMO. 
    396520 
    397521% ------------------------------------------------------------------------------------------------------------- 
    398522%       Linear Bottom Friction 
    399523% ------------------------------------------------------------------------------------------------------------- 
    400 \subsection{Linear Bottom Friction} 
     524\subsection{Linear Bottom Friction (\np{nbotfr} = 1) } 
    401525\label{ZDF_bfr_linear} 
    402526 
    403 %-------------------------------------------nambfr-------------------------------------------------------- 
    404 \begin{flushright} 
    405 (\textbf{namelist} !nbotfr :\textit{ nbotfr = 0, = 1 or = 3}) 
    406 \end{flushright} 
    407  
    408 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): 
     527The linear bottom friction parameterization assumes that the bottom friction  
     528is proportional to the interior velocity (i.e. the velocity of the last model level): 
    409529\begin{equation} \label{Eq_zdfbfr_linear} 
    410 \textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \textbf{U}_h^b 
    411 \end{equation} 
    412 where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom  
    413 ocean layer and $r$ a friction coefficient expressed in m.s$^{-1}$. This  
    414 coefficient is generally estimated by setting a typical decay time $\tau $ in the  
    415 deep ocean, $r = H / \tau$. Commonly accepted values of $\tau$ are of the  
    416 order of 100 to 200 days \citep{Weatherly1984}. A value $\tau^{-1} = 10^{-7}$~s$^{-1}$  
    417 corresponding to 115 days is usually used in quasi-geostrophic models. One may  
    418 consider the linear friction as an approximation of quadratic friction,  
    419 $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}$,  
    420 and assuming an ocean depth $H = 4000$~m, the resulting friction  
    421 coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$. This is the default value used in  
    422 OPA. It corresponds to a decay time scale of 115~days. It can be changed by  
    423 specifying \np{bfric1} (namelist parameter). 
    424  
    425 In the code, the bottom friction is specified by updating the value of the  
    426 vertical eddy coefficient at the bottom level. Indeed, the discrete  
    427 formulation of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that  
    428 $\textbf {U}_h =0$ inside the bottom, leads to 
     530\textbf{F}_h = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b 
     531\end{equation} 
     532where $\textbf{U}_h^b$ is the horizontal velocity vector of the bottom ocean  
     533layer and $r$ is a friction coefficient expressed in m.s$^{-1}$. This coefficient  
     534is generally estimated by setting a typical decay time $\tau$ in the deep ocean,  
     535and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted  
     536values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly1984}.  
     537A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used  
     538in quasi-geostrophic models. One may consider the linear friction as an  
     539approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982},  
     540Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed  
     541of tidal currents of $U_{av} =0.1$~m.s$^{-1}$, and assuming an ocean depth  
     542$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m.s$^{-1}$.  
     543This is the default value used in \NEMO. It corresponds to a decay time scale  
     544of 115~days. It can be changed by specifying \np{bfric1} (namelist parameter). 
     545 
     546In the code, the bottom friction is imposed by updating the value of the  
     547vertical eddy coefficient at the bottom level. Indeed, the discrete formulation  
     548of (\ref{Eq_zdfbfr_linear}) at the last ocean $T-$level, using the fact that  
     549$\textbf {U}_h =0$ below the ocean floor, leads to 
    429550\begin{equation} \label{Eq_zdfbfr_linKz} 
    430551\begin{split} 
    431552A_u^{vm} &= r\;e_{3uw}\\ 
    432 A_v^{vm} &= r\;e_{3uw}\\ 
     553A_v^{vm} &= r\;e_{3vw}\\ 
    433554\end{split} 
    434555\end{equation} 
    435556 
    436 Such an update is done in \mdl{zdfbfr} when \np{nbotfr}=1 and the value of $r$ used is  
    437 \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to set $r=0$ and leads to a free-slip bottom boundary condition,  
    438 while setting \np{nbotfr}=0 imposes $r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the  
    439 background vertical eddy coefficient: a no-slip boundary condition is used.  
     557This update is done in \mdl{zdfbfr} when \np{nbotfr}=1. The value of $r$  
     558used is \np{bfric1}. Setting \np{nbotfr}=3 is equivalent to setting $r=0$ and  
     559leads to a free-slip bottom boundary condition. Setting \np{nbotfr}=0 sets  
     560$r=2\;A_{vb}^{\rm {\bf U}} $, where $A_{vb}^{\rm {\bf U}} $ is the background  
     561vertical eddy coefficient, and a no-slip boundary condition is imposed.  
    440562Note that this latter choice generally leads to an underestimation of the  
    441 bottom friction: for a deepest level thickness of $200~m$ and $A_{vb}^{\rm {\bf U}}  
    442 =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient is only $r=10^{-6}$m.s$^{-1}$. 
    443  
     563bottom friction: for example with a deepest level thickness of $200~m$  
     564and $A_{vb}^{\rm {\bf U}} =10^{-4}$m$^2$.s$^{-1}$, the friction coefficient  
     565is only $r=10^{-6}$m.s$^{-1}$. 
    444566 
    445567% ------------------------------------------------------------------------------------------------------------- 
    446568%       Non-Linear Bottom Friction 
    447569% ------------------------------------------------------------------------------------------------------------- 
    448 \subsection{Non-Linear Bottom Friction} 
     570\subsection{Non-Linear Bottom Friction (\np{nbotfr} = 2)} 
    449571\label{ZDF_bfr_nonlinear} 
    450  
    451 \begin{center} 
    452 (\textbf{namelist} !nbotfr : \textit{nbotfr = 2}) 
    453 \end{center} 
    454572 
    455573The non-linear bottom friction parameterization assumes that the bottom  
    … …  
    460578\end{equation} 
    461579 
    462 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$  
    463 and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been set as  
    464 default value (\np{bfric2} and \np{bfeb2} namelist parameters). 
    465  
    466 As for the linear case, the bottom friction is specified in the code by  
     580with $\textbf{U}_h^b = (u_b\;,\;v_b)$ the horizontal interior velocity ($i.e.$  
     581the horizontal velocity of the bottom ocean layer), $C_D$ a drag coefficient,  
     582and $e_b $ a bottom turbulent kinetic energy due to tides, internal waves  
     583breaking and other short time scale currents. A typical value of the drag  
     584coefficient is $C_D = 10^{-3} $. As an example, the CME experiment  
     585\citep{Treguier1992} uses $C_D = 10^{-3}$ and $e_b = 2.5\;10^{-3}$m$^2$.s$^{-2}$,  
     586while the FRAM experiment \citep{Killworth1992} uses $e_b =0$  
     587and $e_b =2.5\;\;10^{-3}$m$^2$.s$^{-2}$. The FRAM choices have been  
     588set as default values (\np{bfric2} and \np{bfeb2} namelist parameters). 
     589 
     590As for the linear case, the bottom friction is imposed in the code by  
    467591updating the value of the vertical eddy coefficient at the bottom level: 
    468  
    469592\begin{equation} \label{Eq_zdfbfr_nonlinKz} 
    470593\begin{split} 
    471594A_u^{vm} &=C_D\; e_{3uw} \left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^ 
    472595{1/2}\\ 
    473 A_v^{vm} &=C_D\; e_{3uw} \left[  \left(\bar{\bar{u}}^{i+1,j}\right)^2 + v^2 + e_b \right]^{1/2}\\ 
     596A_v^{vm} &=C_D\; e_{3uw} \left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\ 
    474597\end{split} 
    475598\end{equation} 
    476599 
    477600This update is done in \mdl{zdfbfr}. The coefficients that control the strength of the  
    478 non-linear bottom friction are initialized as namelist parameters: ($C_D$= \np{bfri2}, and $e_b$ =\np{bfeb2}). 
    479  
    480 % ================================================================ 
     601non-linear bottom friction are initialized as namelist parameters: $C_D$= \np{bfri2}, and $e_b$ =\np{bfeb2}. 
     602 
     603% ================================================================ 
Note: See TracChangeset for help on using the changeset viewer.