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.
chap_ZDF.tex in branches/2017/dev_merge_2017/DOC/tex_sub – NEMO

source: branches/2017/dev_merge_2017/DOC/tex_sub/chap_ZDF.tex @ 9388

Last change on this file since 9388 was 9388, checked in by nicolasmartin, 6 years ago

Global reorganisation of DOC directory: implementation and configuration of syntax highlighting with 'minted' LaTeX package
New macros have been created to insert different source codes (fortran, xml, c and shell-session) in the document.
Pygments Python package is now required to compile the documentation.

File size: 79.4 KB
Line 
1\documentclass[NEMO_book]{subfiles}
2\begin{document}
3% ================================================================
4% Chapter  Vertical Ocean Physics (ZDF)
5% ================================================================
6\chapter{Vertical Ocean Physics (ZDF)}
7\label{ZDF}
8\minitoc
9
10%gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
11
12
13\newpage
14$\ $\newline    % force a new ligne
15
16
17% ================================================================
18% Vertical Mixing
19% ================================================================
20\section{Vertical Mixing}
21\label{ZDF_zdf}
22
23The discrete form of the ocean subgrid scale physics has been presented in
24\S\ref{TRA_zdf} and \S\ref{DYN_zdf}. At the surface and bottom boundaries,
25the turbulent fluxes of momentum, heat and salt have to be defined. At the
26surface they are prescribed from the surface forcing (see Chap.~\ref{SBC}),
27while at the bottom they are set to zero for heat and salt, unless a geothermal
28flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} 
29defined, see \S\ref{TRA_bbc}), and specified through a bottom friction
30parameterisation for momentum (see \S\ref{ZDF_bfr}).
31
32In this section we briefly discuss the various choices offered to compute
33the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ ,
34$A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$-
35points, respectively (see \S\ref{TRA_zdf} and \S\ref{DYN_zdf}). These
36coefficients can be assumed to be either constant, or a function of the local
37Richardson number, or computed from a turbulent closure model (either
38TKE or GLS formulation). The computation of these coefficients is initialized
39in the \mdl{zdfini} module and performed in the \mdl{zdfric}, \mdl{zdftke} or
40\mdl{zdfgls} modules. The trends due to the vertical momentum and tracer
41diffusion, including the surface forcing, are computed and added to the
42general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.
43These trends can be computed using either a forward time stepping scheme
44(namelist parameter \np{ln\_zdfexp}=true) or a backward time stepping
45scheme (\np{ln\_zdfexp}=false) depending on the magnitude of the mixing
46coefficients, and thus of the formulation used (see \S\ref{STP}).
47
48% -------------------------------------------------------------------------------------------------------------
49%        Constant
50% -------------------------------------------------------------------------------------------------------------
51\subsection{Constant (\protect\key{zdfcst})}
52\label{ZDF_cst}
53%--------------------------------------------namzdf---------------------------------------------------------
54\forfile{../namelists/namzdf}
55%--------------------------------------------------------------------------------------------------------------
56
57Options are defined through the  \ngn{namzdf} namelist variables.
58When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients
59are set to constant values over the whole ocean. This is the crudest way to define
60the vertical ocean physics. It is recommended that this option is only used in
61process studies, not in basin scale simulations. Typical values used in this case are:
62\begin{align*} 
63A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}  \\
64A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
65\end{align*}
66
67These values are set through the \np{rn\_avm0} and \np{rn\_avt0} namelist parameters.
68In all cases, do not use values smaller that those associated with the molecular
69viscosity and diffusivity, that is $\sim10^{-6}~m^2.s^{-1}$ for momentum,
70$\sim10^{-7}~m^2.s^{-1}$ for temperature and $\sim10^{-9}~m^2.s^{-1}$ for salinity.
71
72
73% -------------------------------------------------------------------------------------------------------------
74%        Richardson Number Dependent
75% -------------------------------------------------------------------------------------------------------------
76\subsection{Richardson Number Dependent (\protect\key{zdfric})}
77\label{ZDF_ric}
78
79%--------------------------------------------namric---------------------------------------------------------
80\forfile{../namelists/namzdf_ric}
81%--------------------------------------------------------------------------------------------------------------
82
83When \key{zdfric} is defined, a local Richardson number dependent formulation
84for the vertical momentum and tracer eddy coefficients is set through the  \ngn{namzdf\_ric} 
85namelist variables.The vertical mixing
86coefficients are diagnosed from the large scale variables computed by the model.
87\textit{In situ} measurements have been used to link vertical turbulent activity to
88large scale ocean structures. The hypothesis of a mixing mainly maintained by the
89growth of Kelvin-Helmholtz like instabilities leads to a dependency between the
90vertical eddy coefficients and the local Richardson number ($i.e.$ the
91ratio of stratification to vertical shear). Following \citet{Pacanowski_Philander_JPO81}, the following
92formulation has been implemented:
93\begin{equation} \label{Eq_zdfric}
94   \left\{      \begin{aligned}
95         A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
96         A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
97   \end{aligned}    \right.
98\end{equation}
99where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson
100number, $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),
101$A_b^{vT} $ and $A_b^{vm}$ are the constant background values set as in the
102constant case (see \S\ref{ZDF_cst}), and $A_{ric}^{vT} = 10^{-4}~m^2.s^{-1}$ 
103is the maximum value that can be reached by the coefficient when $Ri\leq 0$,
104$a=5$ and $n=2$. The last three values can be modified by setting the
105\np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively.
106
107A simple mixing-layer model to transfer and dissipate the atmospheric
108 forcings (wind-stress and buoyancy fluxes) can be activated setting
109the \np{ln\_mldw} =.true. in the namelist.
110
111In this case, the local depth of turbulent wind-mixing or "Ekman depth"
112 $h_{e}(x,y,t)$ is evaluated and the vertical eddy coefficients prescribed within this layer.
113
114This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
115\begin{equation}
116         h_{e} = Ek \frac {u^{*}} {f_{0}}    \\
117\end{equation}
118where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis
119parameter.
120
121In this similarity height relationship, the turbulent friction velocity:
122\begin{equation}
123         u^{*} = \sqrt \frac {|\tau|} {\rho_o}     \\
124\end{equation}
125
126is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$.
127The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
128Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
129the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
130
131% -------------------------------------------------------------------------------------------------------------
132%        TKE Turbulent Closure Scheme
133% -------------------------------------------------------------------------------------------------------------
134\subsection{TKE Turbulent Closure Scheme (\protect\key{zdftke})}
135\label{ZDF_tke}
136
137%--------------------------------------------namzdf_tke--------------------------------------------------
138\forfile{../namelists/namzdf_tke}
139%--------------------------------------------------------------------------------------------------------------
140
141The vertical eddy viscosity and diffusivity coefficients are computed from a TKE
142turbulent closure model based on a prognostic equation for $\bar{e}$, the turbulent
143kinetic energy, and a closure assumption for the turbulent length scales. This
144turbulent closure model has been developed by \citet{Bougeault1989} in the
145atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and
146embedded in OPA, the ancestor of NEMO, by \citet{Blanke1993} for equatorial Atlantic
147simulations. Since then, significant modifications have been introduced by
148\citet{Madec1998} in both the implementation and the formulation of the mixing
149length scale. The time evolution of $\bar{e}$ is the result of the production of
150$\bar{e}$ through vertical shear, its destruction through stratification, its vertical
151diffusion, and its dissipation of \citet{Kolmogorov1942} type:
152\begin{equation} \label{Eq_zdftke_e}
153\frac{\partial \bar{e}}{\partial t} =
154\frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
155                    +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
156-K_\rho\,N^2
157+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
158            \;\frac{\partial \bar{e}}{\partial k}} \right]
159- c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
160\end{equation}
161\begin{equation} \label{Eq_zdftke_kz}
162   \begin{split}
163         K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }     \\
164         K_\rho &= A^{vm} / P_{rt}
165   \end{split}
166\end{equation}
167where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2}),
168$l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales,
169$P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity
170and diffusivity coefficients. The constants $C_k =  0.1$ and $C_\epsilon = \sqrt {2} /2$ 
171$\approx 0.7$ are designed to deal with vertical mixing at any depth \citep{Gaspar1990}.
172They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}.
173$P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function
174of the local Richardson number, $R_i$:
175\begin{align*} \label{Eq_prt}
176P_{rt} = \begin{cases}
177                    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}  \\
178                    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}  \\
179                    \ \ 10 &      \text{if $\ 2 \leq R_i$} 
180            \end{cases}
181\end{align*}
182Options are defined through the  \ngn{namzdfy\_tke} namelist variables.
183The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable.
184
185At the sea surface, the value of $\bar{e}$ is prescribed from the wind
186stress field as $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with $e_{bb}$ the \np{rn\_ebb} 
187namelist parameter. The default value of $e_{bb}$ is 3.75. \citep{Gaspar1990}),
188however a much larger value can be used when taking into account the
189surface wave breaking (see below Eq. \eqref{ZDF_Esbc}).
190The bottom value of TKE is assumed to be equal to the value of the level just above.
191The time integration of the $\bar{e}$ equation may formally lead to negative values
192because the numerical scheme does not ensure its positivity. To overcome this
193problem, a cut-off in the minimum value of $\bar{e}$ is used (\np{rn\_emin} 
194namelist parameter). Following \citet{Gaspar1990}, the cut-off value is set
195to $\sqrt{2}/2~10^{-6}~m^2.s^{-2}$. This allows the subsequent formulations
196to match that of \citet{Gargett1984} for the diffusion in the thermocline and
197deep ocean :  $K_\rho = 10^{-3} / N$.
198In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical
199instabilities associated with too weak vertical diffusion. They must be
200specified at least larger than the molecular values, and are set through
201\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
202
203\subsubsection{Turbulent length scale}
204For computational efficiency, the original formulation of the turbulent length
205scales proposed by \citet{Gaspar1990} has been simplified. Four formulations
206are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist
207parameter. The first two are based on the following first order approximation
208\citep{Blanke1993}:
209\begin{equation} \label{Eq_tke_mxl0_1}
210l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
211\end{equation}
212which is valid in a stable stratified region with constant values of the Brunt-
213Vais\"{a}l\"{a} frequency. The resulting length scale is bounded by the distance
214to the surface or to the bottom (\np{nn\_mxl} = 0) or by the local vertical scale factor
215(\np{nn\_mxl} = 1). \citet{Blanke1993} notice that this simplification has two major
216drawbacks: it makes no sense for locally unstable stratification and the
217computation no longer uses all the information contained in the vertical density
218profile. To overcome these drawbacks, \citet{Madec1998} introduces the
219\np{nn\_mxl} = 2 or 3 cases, which add an extra assumption concerning the vertical
220gradient of the computed length scale. So, the length scales are first evaluated
221as in \eqref{Eq_tke_mxl0_1} and then bounded such that:
222\begin{equation} \label{Eq_tke_mxl_constraint}
223\frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
224\qquad \text{with }\  l =  l_k = l_\epsilon
225\end{equation}
226\eqref{Eq_tke_mxl_constraint} means that the vertical variations of the length
227scale cannot be larger than the variations of depth. It provides a better
228approximation of the \citet{Gaspar1990} formulation while being much less
229time consuming. In particular, it allows the length scale to be limited not only
230by the distance to the surface or to the ocean bottom but also by the distance
231to a strongly stratified portion of the water column such as the thermocline
232(Fig.~\ref{Fig_mixing_length}). In order to impose the \eqref{Eq_tke_mxl_constraint} 
233constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$,
234the upward and downward length scales, and evaluate the dissipation and
235mixing length scales as (and note that here we use numerical indexing):
236%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
237\begin{figure}[!t] \begin{center}
238\includegraphics[width=1.00\textwidth]{Fig_mixing_length}
239\caption{ \protect\label{Fig_mixing_length} 
240Illustration of the mixing length computation. }
241\end{center} 
242\end{figure}
243%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
244\begin{equation} \label{Eq_tke_mxl2}
245\begin{aligned}
246 l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
247    \quad &\text{ from $k=1$ to $jpk$ }\ \\
248 l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)   
249    \quad &\text{ from $k=jpk$ to $1$ }\ \\
250\end{aligned}
251\end{equation}
252where $l^{(k)}$ is computed using \eqref{Eq_tke_mxl0_1},
253$i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
254
255In the \np{nn\_mxl}~=~2 case, the dissipation and mixing length scales take the same
256value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the
257\np{nn\_mxl}~=~3 case, the dissipation and mixing turbulent length scales are give
258as in \citet{Gaspar1990}:
259\begin{equation} \label{Eq_tke_mxl_gaspar}
260\begin{aligned}
261& l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }    \\
262& l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)
263\end{aligned}
264\end{equation}
265
266At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist
267parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ 
268where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness
269parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} 
270leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior
271a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ 
272reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
273
274
275\subsubsection{Surface wave breaking parameterization}
276%-----------------------------------------------------------------------%
277Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified
278to include the effect of surface wave breaking energetics. This results in a reduction of summertime
279surface temperature when the mixed layer is relatively shallow. The \citet{Mellor_Blumberg_JPO04} 
280modifications acts on surface length scale and TKE values and air-sea drag coefficient.
281The latter concerns the bulk formulea and is not discussed here.
282
283Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
284\begin{equation}  \label{ZDF_Esbc}
285\bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
286\end{equation}
287where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality
288which depends on the ''wave age'', ranging from 57 for mature waves to 146 for
289younger waves \citep{Mellor_Blumberg_JPO04}.
290The boundary condition on the turbulent length scale follows the Charnock's relation:
291\begin{equation} \label{ZDF_Lsbc}
292l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
293\end{equation}
294where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
295\citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by \citet{Stacey_JPO99}
296citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value.
297As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,
298with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds
299to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} 
300as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value.
301Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters)
302is applied on surface $\bar{e}$ value.
303
304
305\subsubsection{Langmuir cells}
306%--------------------------------------%
307Langmuir circulations (LC) can be described as ordered large-scale vertical motions
308in the surface layer of the oceans. Although LC have nothing to do with convection,
309the circulation pattern is rather similar to so-called convective rolls in the atmospheric
310boundary layer. The detailed physics behind LC is described in, for example,
311\citet{Craik_Leibovich_JFM76}. The prevailing explanation is that LC arise from
312a nonlinear interaction between the Stokes drift and wind drift currents.
313
314Here we introduced in the TKE turbulent closure the simple parameterization of
315Langmuir circulations proposed by \citep{Axell_JGR02} for a $k-\epsilon$ turbulent closure.
316The parameterization, tuned against large-eddy simulation, includes the whole effect
317of LC in an extra source terms of TKE, $P_{LC}$.
318The presence of $P_{LC}$ in \eqref{Eq_zdftke_e}, the TKE equation, is controlled
319by setting \np{ln\_lc} to \textit{true} in the namtke namelist.
320 
321By making an analogy with the characteristic convective velocity scale
322($e.g.$, \citet{D'Alessio_al_JPO98}), $P_{LC}$ is assumed to be :
323\begin{equation}
324P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}
325\end{equation}
326where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
327With no information about the wave field, $w_{LC}$ is assumed to be proportional to
328the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
329\footnote{Following \citet{Li_Garrett_JMR93}, the surface Stoke drift velocity
330may be expressed as $u_s =  0.016 \,|U_{10m}|$. Assuming an air density of
331$\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of $1.5~10^{-3}$ give the expression
332used of $u_s$ as a function of the module of surface stress}.
333For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as
334at a finite depth $H_{LC}$ (which is often close to the mixed layer depth), and simply
335varies as a sine function in between (a first-order profile for the Langmuir cell structures).
336The resulting expression for $w_{LC}$ is :
337\begin{equation}
338w_{LC}  = \begin{cases}
339                   c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$}    \\
340                   0                             &      \text{otherwise} 
341                 \end{cases}
342\end{equation}
343where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise
344to fit LES data. The chosen value yields maximum vertical velocities $w_{LC}$ of the order
345of a few centimeters per second. The value of $c_{LC}$ is set through the \np{rn\_lc} 
346namelist parameter, having in mind that it should stay between 0.15 and 0.54 \citep{Axell_JGR02}.
347
348The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
349$H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift
350can reach on its own by converting its kinetic energy to potential energy, according to
351\begin{equation}
352- \int_{-H_{LC}}^0 { N^2\;\;dz} = \frac{1}{2} u_s^2
353\end{equation}
354
355
356\subsubsection{Mixing just below the mixed layer}
357%--------------------------------------------------------------%
358
359Vertical mixing parameterizations commonly used in ocean general circulation models
360tend to produce mixed-layer depths that are too shallow during summer months and windy conditions.
361This bias is particularly acute over the Southern Ocean.
362To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme  \cite{Rodgers_2014}.
363The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations,
364but rather is meant to account for observed processes that affect the density structure of
365the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme
366($i.e.$ near-inertial oscillations and ocean swells and waves).
367
368When using this parameterization ($i.e.$ when \np{nn\_etau}~=~1), the TKE input to the ocean ($S$)
369imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized
370by \eqref{ZDF_Esbc} the standard TKE surface boundary condition, plus a depth depend one given by:
371\begin{equation}  \label{ZDF_Ehtau}
372S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} 
373\end{equation}
374where
375$z$ is the depth, 
376$e_s$ is TKE surface boundary condition,
377$f_r$ is the fraction of the surface TKE that penetrate in the ocean,
378$h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration,
379and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely
380covered by sea-ice).
381The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter.
382The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}~=~0)
383or a latitude dependent value (varying from 0.5~m at the Equator to a maximum value of 30~m
384at high latitudes (\np{nn\_etau}~=~1).
385
386Note that two other option existe, \np{nn\_etau}~=~2, or 3. They correspond to applying
387\eqref{ZDF_Ehtau} only at the base of the mixed layer, or to using the high frequency part
388of the stress to evaluate the fraction of TKE that penetrate the ocean.
389Those two options are obsolescent features introduced for test purposes.
390They will be removed in the next release.
391
392
393
394% from Burchard et al OM 2008 :
395% the most critical process not reproduced by statistical turbulence models is the activity of
396% internal waves and their interaction with turbulence. After the Reynolds decomposition,
397% internal waves are in principle included in the RANS equations, but later partially
398% excluded by the hydrostatic assumption and the model resolution.
399% Thus far, the representation of internal wave mixing in ocean models has been relatively crude
400% (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
401
402
403
404% -------------------------------------------------------------------------------------------------------------
405%        TKE discretization considerations
406% -------------------------------------------------------------------------------------------------------------
407\subsection{TKE discretization considerations (\protect\key{zdftke})}
408\label{ZDF_tke_ene}
409
410%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
411\begin{figure}[!t]   \begin{center}
412\includegraphics[width=1.00\textwidth]{Fig_ZDF_TKE_time_scheme}
413\caption{ \protect\label{Fig_TKE_time_scheme} 
414Illustration of the TKE time integration and its links to the momentum and tracer time integration. }
415\end{center} 
416\end{figure}
417%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
418
419The production of turbulence by vertical shear (the first term of the right hand side
420of \eqref{Eq_zdftke_e}) should balance the loss of kinetic energy associated with
421the vertical momentum diffusion (first line in \eqref{Eq_PE_zdf}). To do so a special care
422have to be taken for both the time and space discretization of the TKE equation
423\citep{Burchard_OM02,Marsaleix_al_OM08}.
424
425Let us first address the time stepping issue. Fig.~\ref{Fig_TKE_time_scheme} shows
426how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays
427with the one-level forward time stepping of TKE equation. With this framework, the total loss
428of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is
429obtained by multiplying this quantity by $u^t$ and summing the result vertically:   
430\begin{equation} \label{Eq_energ1}
431\begin{split}
432\int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
433&= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}         
434 - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
435\end{split}
436\end{equation}
437Here, the vertical diffusion of momentum is discretized backward in time
438with a coefficient, $K_m$, known at time $t$ (Fig.~\ref{Fig_TKE_time_scheme}),
439as it is required when using the TKE scheme (see \S\ref{STP_forward_imp}).
440The first term of the right hand side of \eqref{Eq_energ1} represents the kinetic energy
441transfer at the surface (atmospheric forcing) and at the bottom (friction effect).
442The second term is always negative. It is the dissipation rate of kinetic energy,
443and thus minus the shear production rate of $\bar{e}$. \eqref{Eq_energ1} 
444implies that, to be energetically consistent, the production rate of $\bar{e}$ 
445used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
446${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$ (and not by the more straightforward
447$K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
448
449A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
450(second term of the right hand side of \eqref{Eq_zdftke_e}). This term
451must balance the input of potential energy resulting from vertical mixing.
452The rate of change of potential energy (in 1D for the demonstration) due vertical
453mixing is obtained by multiplying vertical density diffusion
454tendency by $g\,z$ and and summing the result vertically:
455\begin{equation} \label{Eq_energ2}
456\begin{split}
457\int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
458&= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta} 
459   - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
460&= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
461+ \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
462\end{split}
463\end{equation}
464where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
465The first term of the right hand side of \eqref{Eq_energ2}  is always zero
466because there is no diffusive flux through the ocean surface and bottom).
467The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
468Therefore \eqref{Eq_energ1} implies that, to be energetically consistent, the product
469${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \eqref{Eq_zdftke_e}, the TKE equation.
470
471Let us now address the space discretization issue.
472The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity
473components are in the centre of the side faces of a $t$-box in staggered C-grid
474(Fig.\ref{Fig_cell}). A space averaging is thus required to obtain the shear TKE production term.
475By redoing the \eqref{Eq_energ1} in the 3D case, it can be shown that the product of
476eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
477Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into
478account.
479
480The above energetic considerations leads to
481the following final discrete form for the TKE equation:
482\begin{equation} \label{Eq_zdftke_ene}
483\begin{split}
484\frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv 
485\Biggl\{ \Biggr.
486  &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} } 
487                                                                              \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
488+&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} } 
489                                                                               \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j} 
490\Biggr. \Biggr\}   \\
491%
492- &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
493%
494+&\frac{1}{{e_3w}^{t+\rdt}}  \;\delta_{k+1/2} \left[   {K_m}^{t-\rdt} \,\frac{\delta_{k}[(\bar{e})^{t+\rdt}]} {{e_3w}^{t+\rdt}}   \right]   \\
495%
496- &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
497\end{split}
498\end{equation}
499where the last two terms in \eqref{Eq_zdftke_ene} (vertical diffusion and Kolmogorov dissipation)
500are time stepped using a backward scheme (see\S\ref{STP_forward_imp}).
501Note that the Kolmogorov term has been linearized in time in order to render
502the implicit computation possible. The restart of the TKE scheme
503requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in
504the right hand side of \eqref{Eq_zdftke_ene}. For the latter, it is in fact
505the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
506
507% -------------------------------------------------------------------------------------------------------------
508%        GLS Generic Length Scale Scheme
509% -------------------------------------------------------------------------------------------------------------
510\subsection{GLS Generic Length Scale (\protect\key{zdfgls})}
511\label{ZDF_gls}
512
513%--------------------------------------------namzdf_gls---------------------------------------------------------
514\forfile{../namelists/namzdf_gls}
515%--------------------------------------------------------------------------------------------------------------
516
517The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on
518two prognostic equations: one for the turbulent kinetic energy $\bar {e}$, and another
519for the generic length scale, $\psi$ \citep{Umlauf_Burchard_JMS03, Umlauf_Burchard_CSR05}.
520This later variable is defined as : $\psi = {C_{0\mu}}^{p} \ {\bar{e}}^{m} \ l^{n}$,
521where the triplet $(p, m, n)$ value given in Tab.\ref{Tab_GLS} allows to recover
522a number of well-known turbulent closures ($k$-$kl$ \citep{Mellor_Yamada_1982},
523$k$-$\epsilon$ \citep{Rodi_1987}, $k$-$\omega$ \citep{Wilcox_1988} 
524among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}).
525The GLS scheme is given by the following set of equations:
526\begin{equation} \label{Eq_zdfgls_e}
527\frac{\partial \bar{e}}{\partial t} =
528\frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
529                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
530-K_\rho \,N^2
531+\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
532- \epsilon
533\end{equation}
534
535\begin{equation} \label{Eq_zdfgls_psi}
536   \begin{split}
537\frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
538\frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
539                                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
540- C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
541&+\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
542                                  \;\frac{\partial \psi}{\partial k}} \right]\;
543   \end{split}
544\end{equation}
545
546\begin{equation} \label{Eq_zdfgls_kz}
547   \begin{split}
548         K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
549         K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
550   \end{split}
551\end{equation}
552
553\begin{equation} \label{Eq_zdfgls_eps}
554{\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
555\end{equation}
556where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \S\ref{TRA_bn2})
557and $\epsilon$ the dissipation rate.
558The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$)
559depends of the choice of the turbulence model. Four different turbulent models are pre-defined
560(Tab.\ref{Tab_GLS}). They are made available through the \np{nn\_clo} namelist parameter.
561
562%--------------------------------------------------TABLE--------------------------------------------------
563\begin{table}[htbp]  \begin{center}
564%\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
565\begin{tabular}{ccccc}
566                         &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\ 
567%                        & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\ 
568\hline  \hline 
569\np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\ 
570\hline 
571$( p , n , m )$          &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
572$\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
573$\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
574$C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
575$C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
576$C_3$              &      1.           &     1.              &      1.                &       1.           \\
577$F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
578\hline
579\hline
580\end{tabular}
581\caption{   \protect\label{Tab_GLS} 
582Set of predefined GLS parameters, or equivalently predefined turbulence models available
583with \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls} .}
584\end{center}   \end{table}
585%--------------------------------------------------------------------------------------------------------------
586
587In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force
588the convergence of the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length)
589value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ 
590are calculated from stability function proposed by \citet{Galperin_al_JAS88}, or by \citet{Kantha_Clayson_1994} 
591or one of the two functions suggested by \citet{Canuto_2001}  (\np{nn\_stab\_func} = 0, 1, 2 or 3, resp.).
592The value of $C_{0\mu}$ depends of the choice of the stability function.
593
594The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated
595thanks to Dirichlet or Neumann condition through \np{nn\_tkebc\_surf} and \np{nn\_tkebc\_bot}, resp.
596As for TKE closure , the wave effect on the mixing is considered when \np{ln\_crban}~=~true
597\citep{Craig_Banner_JPO94, Mellor_Blumberg_JPO04}. The \np{rn\_crban} namelist parameter
598is $\alpha_{CB}$ in \eqref{ZDF_Esbc} and \np{rn\_charn} provides the value of $\beta$ in \eqref{ZDF_Lsbc}.
599
600The $\psi$ equation is known to fail in stably stratified flows, and for this reason
601almost all authors apply a clipping of the length scale as an \textit{ad hoc} remedy.
602With this clipping, the maximum permissible length scale is determined by
603$l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. A value of $c_{lim} = 0.53$ is often used
604\citep{Galperin_al_JAS88}. \cite{Umlauf_Burchard_CSR05} show that the value of
605the clipping factor is of crucial importance for the entrainment depth predicted in
606stably stratified situations, and that its value has to be chosen in accordance
607with the algebraic model for the turbulent fluxes. The clipping is only activated
608if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
609
610The time and space discretization of the GLS equations follows the same energetic
611consideration as for the TKE case described in \S\ref{ZDF_tke_ene}  \citep{Burchard_OM02}.
612Examples of performance of the 4 turbulent closure scheme can be found in \citet{Warner_al_OM05}.
613
614% -------------------------------------------------------------------------------------------------------------
615%        OSM OSMOSIS BL Scheme
616% -------------------------------------------------------------------------------------------------------------
617\subsection{OSM OSMOSIS Boundary Layer scheme (\protect\key{zdfosm})}
618\label{ZDF_osm}
619
620%--------------------------------------------namzdf_osm---------------------------------------------------------
621\forfile{../namelists/namzdf_osm}
622%--------------------------------------------------------------------------------------------------------------
623
624The OSMOSIS turbulent closure scheme is based on......   TBC
625
626% ================================================================
627% Convection
628% ================================================================
629\section{Convection}
630\label{ZDF_conv}
631
632%--------------------------------------------namzdf--------------------------------------------------------
633\forfile{../namelists/namzdf}
634%--------------------------------------------------------------------------------------------------------------
635
636Static instabilities (i.e. light potential densities under heavy ones) may
637occur at particular ocean grid points. In nature, convective processes
638quickly re-establish the static stability of the water column. These
639processes have been removed from the model via the hydrostatic
640assumption so they must be parameterized. Three parameterisations
641are available to deal with convective processes: a non-penetrative
642convective adjustment or an enhanced vertical diffusion, or/and the
643use of a turbulent closure scheme.
644
645% -------------------------------------------------------------------------------------------------------------
646%       Non-Penetrative Convective Adjustment
647% -------------------------------------------------------------------------------------------------------------
648\subsection   [Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}) ]
649         {Non-Penetrative Convective Adjustment (\protect\np{ln\_tranpc}=.true.) }
650\label{ZDF_npc}
651
652%--------------------------------------------namzdf--------------------------------------------------------
653\forfile{../namelists/namzdf}
654%--------------------------------------------------------------------------------------------------------------
655
656%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
657\begin{figure}[!htb]    \begin{center}
658\includegraphics[width=0.90\textwidth]{Fig_npc}
659\caption{  \protect\label{Fig_npc} 
660Example of an unstable density profile treated by the non penetrative
661convective adjustment algorithm. $1^{st}$ step: the initial profile is checked from
662the surface to the bottom. It is found to be unstable between levels 3 and 4.
663They are mixed. The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5
664are mixed. The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are
665mixed. The $1^{st}$ step ends since the density profile is then stable below
666the level 3. $2^{nd}$ step: the new $\rho$ profile is checked following the same
667procedure as in $1^{st}$ step: levels 2 to 5 are mixed. The new density profile
668is checked. It is found stable: end of algorithm.}
669\end{center}   \end{figure}
670%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
671
672Options are defined through the  \ngn{namzdf} namelist variables.
673The non-penetrative convective adjustment is used when \np{ln\_zdfnpc}~=~\textit{true}.
674It is applied at each \np{nn\_npc} time step and mixes downwards instantaneously
675the statically unstable portion of the water column, but only until the density
676structure becomes neutrally stable ($i.e.$ until the mixed portion of the water
677column has \textit{exactly} the density of the water just below) \citep{Madec_al_JPO91}.
678The associated algorithm is an iterative process used in the following way
679(Fig. \ref{Fig_npc}): starting from the top of the ocean, the first instability is
680found. Assume in the following that the instability is located between levels
681$k$ and $k+1$. The temperature and salinity in the two levels are
682vertically mixed, conserving the heat and salt contents of the water column.
683The new density is then computed by a linear approximation. If the new
684density profile is still unstable between levels $k+1$ and $k+2$, levels $k$,
685$k+1$ and $k+2$ are then mixed. This process is repeated until stability is
686established below the level $k$ (the mixing process can go down to the
687ocean bottom). The algorithm is repeated to check if the density profile
688between level $k-1$ and $k$ is unstable and/or if there is no deeper instability.
689
690This algorithm is significantly different from mixing statically unstable levels
691two by two. The latter procedure cannot converge with a finite number
692of iterations for some vertical profiles while the algorithm used in \NEMO 
693converges for any profile in a number of iterations which is less than the
694number of vertical levels. This property is of paramount importance as
695pointed out by \citet{Killworth1989}: it avoids the existence of permanent
696and unrealistic static instabilities at the sea surface. This non-penetrative
697convective algorithm has been proved successful in studies of the deep
698water formation in the north-western Mediterranean Sea
699\citep{Madec_al_JPO91, Madec_al_DAO91, Madec_Crepon_Bk91}.
700
701The current implementation has been modified in order to deal with any non linear
702equation of seawater (L. Brodeau, personnal communication).
703Two main differences have been introduced compared to the original algorithm:
704$(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency
705(not the the difference in potential density) ;
706$(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients
707are vertically mixed in the same way their temperature and salinity has been mixed.
708These two modifications allow the algorithm to perform properly and accurately
709with TEOS10 or EOS-80 without having to recompute the expansion coefficients at each
710mixing iteration.
711
712% -------------------------------------------------------------------------------------------------------------
713%       Enhanced Vertical Diffusion
714% -------------------------------------------------------------------------------------------------------------
715\subsection   [Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd})]
716              {Enhanced Vertical Diffusion (\protect\np{ln\_zdfevd}=true)}
717\label{ZDF_evd}
718
719%--------------------------------------------namzdf--------------------------------------------------------
720\forfile{../namelists/namzdf}
721%--------------------------------------------------------------------------------------------------------------
722
723Options are defined through the  \ngn{namzdf} namelist variables.
724The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true.
725In this case, the vertical eddy mixing coefficients are assigned very large values
726(a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable
727($i.e.$ when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative)
728\citep{Lazar_PhD97, Lazar_al_JPO99}. This is done either on tracers only
729(\np{nn\_evdm}=0) or on both momentum and tracers (\np{nn\_evdm}=1).
730
731In practice, where $N^2\leq 10^{-12}$, $A_T^{vT}$ and $A_T^{vS}$, and
732if \np{nn\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
733values also, are set equal to the namelist parameter \np{rn\_avevd}. A typical value
734for $rn\_avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of
735convective processes is less time consuming than the convective adjustment
736algorithm presented above when mixing both tracers and momentum in the
737case of static instabilities. It requires the use of an implicit time stepping on
738vertical diffusion terms (i.e. \np{ln\_zdfexp}=false).
739
740Note that the stability test is performed on both \textit{before} and \textit{now} 
741values of $N^2$. This removes a potential source of divergence of odd and
742even time step in a leapfrog environment \citep{Leclair_PhD2010} (see \S\ref{STP_mLF}).
743
744% -------------------------------------------------------------------------------------------------------------
745%       Turbulent Closure Scheme
746% -------------------------------------------------------------------------------------------------------------
747\subsection[Turbulent Closure Scheme (\protect\key{zdf\{tke, gls, osm\}})]{Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})}
748\label{ZDF_tcs}
749
750The turbulent closure scheme presented in \S\ref{ZDF_tke} and \S\ref{ZDF_gls} 
751(\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically
752unstable density profiles. In such a case, the term corresponding to the
753destruction of turbulent kinetic energy through stratification in \eqref{Eq_zdftke_e} 
754or \eqref{Eq_zdfgls_e} becomes a source term, since $N^2$ is negative.
755It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring
756$A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values
757restore the static stability of the water column in a way similar to that of the
758enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However,
759in the vicinity of the sea surface (first ocean layer), the eddy coefficients
760computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$,
761because the mixing length scale is bounded by the distance to the sea surface.
762It can thus be useful to combine the enhanced vertical
763diffusion with the turbulent closure scheme, $i.e.$ setting the \np{ln\_zdfnpc} 
764namelist parameter to true and defining the turbulent closure CPP key all together.
765
766The KPP turbulent closure scheme already includes enhanced vertical diffusion
767in the case of convection, as governed by the variables $bvsqcon$ and $difcon$ 
768found in \mdl{zdfkpp}, therefore \np{ln\_zdfevd}=false should be used with the KPP
769scheme. %gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
770
771% ================================================================
772% Double Diffusion Mixing
773% ================================================================
774\section  [Double Diffusion Mixing (\protect\key{zdfddm})]
775      {Double Diffusion Mixing (\protect\key{zdfddm})}
776\label{ZDF_ddm}
777
778%-------------------------------------------namzdf_ddm-------------------------------------------------
779%\forfile{../namelists/namzdf_ddm}
780%--------------------------------------------------------------------------------------------------------------
781
782Options are defined through the  \ngn{namzdf\_ddm} namelist variables.
783Double diffusion occurs when relatively warm, salty water overlies cooler, fresher
784water, or vice versa. The former condition leads to salt fingering and the latter
785to diffusive convection. Double-diffusive phenomena contribute to diapycnal
786mixing in extensive regions of the ocean.  \citet{Merryfield1999} include a
787parameterisation of such phenomena in a global ocean model and show that
788it leads to relatively minor changes in circulation but exerts significant regional
789influences on temperature and salinity. This parameterisation has been
790introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key.
791
792Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
793\begin{align*} % \label{Eq_zdfddm_Kz}
794    &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT}  \\
795    &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
796\end{align*}
797where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
798and $o$ by processes other than double diffusion. The rates of double-diffusive
799mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$,
800where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline
801contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt
802fingering, we adopt the diapycnal diffusivities suggested by Schmitt (1981):
803\begin{align} \label{Eq_zdfddm_f}
804A_f^{vS} &=    \begin{cases}
805   \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
806   0                              &\text{otherwise} 
807            \end{cases}   
808\\           \label{Eq_zdfddm_f_T}
809A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
810\end{align}
811
812%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
813\begin{figure}[!t]   \begin{center}
814\includegraphics[width=0.99\textwidth]{Fig_zdfddm}
815\caption{  \protect\label{Fig_zdfddm}
816From \citet{Merryfield1999} : (a) Diapycnal diffusivities $A_f^{vT}$ 
817and $A_f^{vS}$ for temperature and salt in regions of salt fingering. Heavy
818curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves
819$A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and
820$A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy
821curves denote the Federov parameterisation and thin curves the Kelley
822parameterisation. The latter is not implemented in \NEMO. }
823\end{center}    \end{figure}
824%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
825
826The factor 0.7 in \eqref{Eq_zdfddm_f_T} reflects the measured ratio
827$\alpha F_T /\beta F_S \approx  0.7$ of buoyancy flux of heat to buoyancy
828flux of salt ($e.g.$, \citet{McDougall_Taylor_JMR84}). Following  \citet{Merryfield1999},
829we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
830
831To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used:
832\begin{align}  \label{Eq_zdfddm_d}
833A_d^{vT} &=    \begin{cases}
834   1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
835                           &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
836   0                       &\text{otherwise} 
837            \end{cases}   
838\\          \label{Eq_zdfddm_d_S}
839A_d^{vS} &=    \begin{cases}
840   A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right)
841                           &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
842   A_d^{vT} \ 0.15 \ R_\rho
843                           &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
844   0                       &\text{otherwise} 
845            \end{cases}   
846\end{align}
847
848The dependencies of \eqref{Eq_zdfddm_f} to \eqref{Eq_zdfddm_d_S} on $R_\rho$ 
849are illustrated in Fig.~\ref{Fig_zdfddm}. Implementing this requires computing
850$R_\rho$ at each grid point on every time step. This is done in \mdl{eosbn2} at the
851same time as $N^2$ is computed. This avoids duplication in the computation of
852$\alpha$ and $\beta$ (which is usually quite expensive).
853
854% ================================================================
855% Bottom Friction
856% ================================================================
857\section  [Bottom and Top Friction (\textit{zdfbfr})]   {Bottom and Top Friction (\protect\mdl{zdfbfr} module)}
858\label{ZDF_bfr}
859
860%--------------------------------------------nambfr--------------------------------------------------------
861%\forfile{../namelists/nambfr}
862%--------------------------------------------------------------------------------------------------------------
863
864Options to define the top and bottom friction are defined through the  \ngn{nambfr} namelist variables.
865The bottom friction represents the friction generated by the bathymetry.
866The top friction represents the friction generated by the ice shelf/ocean interface.
867As the friction processes at the top and bottom are treated in similar way,
868only the bottom friction is described in detail below.
869
870
871Both the surface momentum flux (wind stress) and the bottom momentum
872flux (bottom friction) enter the equations as a condition on the vertical
873diffusive flux. For the bottom boundary layer, one has:
874\begin{equation} \label{Eq_zdfbfr_flux}
875A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
876\end{equation}
877where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum
878outside the logarithmic turbulent boundary layer (thickness of the order of
8791~m in the ocean). How ${\cal F}_h^{\textbf U}$ influences the interior depends on the
880vertical resolution of the model near the bottom relative to the Ekman layer
881depth. For example, in order to obtain an Ekman layer depth
882$d = \sqrt{2\;A^{vm}} / f = 50$~m, one needs a vertical diffusion coefficient
883$A^{vm} = 0.125$~m$^2$s$^{-1}$ (for a Coriolis frequency
884$f = 10^{-4}$~m$^2$s$^{-1}$). With a background diffusion coefficient
885$A^{vm} = 10^{-4}$~m$^2$s$^{-1}$, the Ekman layer depth is only 1.4~m.
886When the vertical mixing coefficient is this small, using a flux condition is
887equivalent to entering the viscous forces (either wind stress or bottom friction)
888as a body force over the depth of the top or bottom model layer. To illustrate
889this, consider the equation for $u$ at $k$, the last ocean level:
890\begin{equation} \label{Eq_zdfbfr_flux2}
891\frac{\partial u_k}{\partial t} = \frac{1}{e_{3u}} \left[ \frac{A_{uw}^{vm}}{e_{3uw}} \delta_{k+1/2}\;[u] - {\cal F}^u_h \right] \approx - \frac{{\cal F}^u_{h}}{e_{3u}}
892\end{equation}
893If the bottom layer thickness is 200~m, the Ekman transport will
894be distributed over that depth. On the other hand, if the vertical resolution
895is high (1~m or less) and a turbulent closure model is used, the turbulent
896Ekman layer will be represented explicitly by the model. However, the
897logarithmic layer is never represented in current primitive equation model
898applications: it is \emph{necessary} to parameterize the flux ${\cal F}^u_h $.
899Two choices are available in \NEMO: a linear and a quadratic bottom friction.
900Note that in both cases, the rotation between the interior velocity and the
901bottom friction is neglected in the present release of \NEMO.
902
903In the code, the bottom friction is imposed by adding the trend due to the bottom
904friction to the general momentum trend in \mdl{dynbfr}. For the time-split surface
905pressure gradient algorithm, the momentum trend due to the barotropic component
906needs to be handled separately. For this purpose it is convenient to compute and
907store coefficients which can be simply combined with bottom velocities and geometric
908values to provide the momentum trend due to bottom friction.
909These coefficients are computed in \mdl{zdfbfr} and generally take the form
910$c_b^{\textbf U}$ where:
911\begin{equation} \label{Eq_zdfbfr_bdef}
912\frac{\partial {\textbf U_h}}{\partial t} =
913  - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
914\end{equation}
915where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
916
917% -------------------------------------------------------------------------------------------------------------
918%       Linear Bottom Friction
919% -------------------------------------------------------------------------------------------------------------
920\subsection{Linear Bottom Friction (\protect\np{nn\_botfr} = 0 or 1) }
921\label{ZDF_bfr_linear}
922
923The linear bottom friction parameterisation (including the special case
924of a free-slip condition) assumes that the bottom friction
925is proportional to the interior velocity (i.e. the velocity of the last
926model level):
927\begin{equation} \label{Eq_zdfbfr_linear}
928{\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
929\end{equation}
930where $r$ is a friction coefficient expressed in ms$^{-1}$.
931This coefficient is generally estimated by setting a typical decay time
932$\tau$ in the deep ocean,
933and setting $r = H / \tau$, where $H$ is the ocean depth. Commonly accepted
934values of $\tau$ are of the order of 100 to 200 days \citep{Weatherly_JMR84}.
935A value $\tau^{-1} = 10^{-7}$~s$^{-1}$ equivalent to 115 days, is usually used
936in quasi-geostrophic models. One may consider the linear friction as an
937approximation of quadratic friction, $r \approx 2\;C_D\;U_{av}$ (\citet{Gill1982},
938Eq. 9.6.6). For example, with a drag coefficient $C_D = 0.002$, a typical speed
939of tidal currents of $U_{av} =0.1$~m\;s$^{-1}$, and assuming an ocean depth
940$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
941This is the default value used in \NEMO. It corresponds to a decay time scale
942of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
943
944For the linear friction case the coefficients defined in the general
945expression \eqref{Eq_zdfbfr_bdef} are:
946\begin{equation} \label{Eq_zdfbfr_linbfr_b}
947\begin{split}
948 c_b^u &= - r\\
949 c_b^v &= - r\\
950\end{split}
951\end{equation}
952When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}.
953Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip
954bottom boundary condition. These values are assigned in \mdl{zdfbfr}.
955From v3.2 onwards there is support for local enhancement of these values
956via an externally defined 2D mask array (\np{ln\_bfr2d}=true) given
957in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1.
958Locations with a non-zero mask value will have the friction coefficient increased
959by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
960
961% -------------------------------------------------------------------------------------------------------------
962%       Non-Linear Bottom Friction
963% -------------------------------------------------------------------------------------------------------------
964\subsection{Non-Linear Bottom Friction (\protect\np{nn\_botfr} = 2)}
965\label{ZDF_bfr_nonlinear}
966
967The non-linear bottom friction parameterisation assumes that the bottom
968friction is quadratic:
969\begin{equation} \label{Eq_zdfbfr_nonlinear}
970{\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h
971}{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b
972\end{equation}
973where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy
974due to tides, internal waves breaking and other short time scale currents.
975A typical value of the drag coefficient is $C_D = 10^{-3} $. As an example,
976the CME experiment \citep{Treguier_JGR92} uses $C_D = 10^{-3}$ and
977$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} 
978uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
979The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} 
980namelist parameters).
981
982As for the linear case, the bottom friction is imposed in the code by
983adding the trend due to the bottom friction to the general momentum trend
984in \mdl{dynbfr}.
985For the non-linear friction case the terms
986computed in \mdl{zdfbfr}  are:
987\begin{equation} \label{Eq_zdfbfr_nonlinbfr}
988\begin{split}
989 c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\
990 c_b^v &= - \; C_D\;\left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
991\end{split}
992\end{equation}
993
994The coefficients that control the strength of the non-linear bottom friction are
995initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.
996Note for applications which treat tides explicitly a low or even zero value of
997\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible
998via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  This works in the same way
999as for the linear bottom friction case with non-zero masked locations increased by
1000$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
1001
1002% -------------------------------------------------------------------------------------------------------------
1003%       Bottom Friction Log-layer
1004% -------------------------------------------------------------------------------------------------------------
1005\subsection[Log-layer Bottom Friction enhancement (\protect\np{ln\_loglayer} = .true.)]{Log-layer Bottom Friction enhancement (\protect\np{nn\_botfr} = 2, \protect\np{ln\_loglayer} = .true.)}
1006\label{ZDF_bfr_loglayer}
1007
1008In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally
1009enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no
1010longer constant but is related to the thickness of the last wet layer in each column by:
1011
1012\begin{equation}
1013C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
1014\end{equation}
1015
1016\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness
1017length provided via the namelist.
1018
1019For stability, the drag coefficient is bounded such that it is kept greater or equal to
1020the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional
1021namelist parameter: \np{rn\_bfri2\_max}, i.e.:
1022
1023\begin{equation}
1024rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
1025\end{equation}
1026
1027\noindent Note also that a log-layer enhancement can also be applied to the top boundary
1028friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the
1029relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2}
1030and \np{rn\_tfri2\_max}.
1031
1032% -------------------------------------------------------------------------------------------------------------
1033%       Bottom Friction stability
1034% -------------------------------------------------------------------------------------------------------------
1035\subsection{Bottom Friction stability considerations}
1036\label{ZDF_bfr_stability}
1037
1038Some care needs to exercised over the choice of parameters to ensure that the
1039implementation of bottom friction does not induce numerical instability. For
1040the purposes of stability analysis, an approximation to \eqref{Eq_zdfbfr_flux2}
1041is:
1042\begin{equation} \label{Eqn_bfrstab}
1043\begin{split}
1044 \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
1045               &= -\frac{ru}{e_{3u}}\;2\rdt\\
1046\end{split}
1047\end{equation}
1048\noindent where linear bottom friction and a leapfrog timestep have been assumed.
1049To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have:
1050\begin{equation}
1051 |\Delta u| < \;|u|
1052\end{equation}
1053\noindent which, using \eqref{Eqn_bfrstab}, gives:
1054\begin{equation}
1055r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
1056\end{equation}
1057This same inequality can also be derived in the non-linear bottom friction case
1058if a velocity of 1 m.s$^{-1}$ is assumed. Alternatively, this criterion can be
1059rearranged to suggest a minimum bottom box thickness to ensure stability:
1060\begin{equation}
1061e_{3u} > 2\;r\;\rdt
1062\end{equation}
1063\noindent which it may be necessary to impose if partial steps are being used.
1064For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then
1065$e_{3u}$ should be greater than 3.6 m. For most applications, with physically
1066sensible parameters these restrictions should not be of concern. But
1067caution may be necessary if attempts are made to locally enhance the bottom
1068friction parameters.
1069To ensure stability limits are imposed on the bottom friction coefficients both during
1070initialisation and at each time step. Checks at initialisation are made in \mdl{zdfbfr} 
1071(assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
1072The number of breaches of the stability criterion are reported as well as the minimum
1073and maximum values that have been set. The criterion is also checked at each time step,
1074using the actual velocity, in \mdl{dynbfr}. Values of the bottom friction coefficient are
1075reduced as necessary to ensure stability; these changes are not reported.
1076
1077Limits on the bottom friction coefficient are not imposed if the user has elected to
1078handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential
1079breaches of the explicit stability criterion are still reported for information purposes.
1080
1081% -------------------------------------------------------------------------------------------------------------
1082%       Implicit Bottom Friction
1083% -------------------------------------------------------------------------------------------------------------
1084\subsection[Implicit Bottom Friction (\protect\np{ln\_bfrimp})]{Implicit Bottom Friction (\protect\np{ln\_bfrimp}$=$\textit{T})}
1085\label{ZDF_bfr_imp}
1086
1087An optional implicit form of bottom friction has been implemented to improve
1088model stability. We recommend this option for shelf sea and coastal ocean applications, especially
1089for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} 
1090to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false} 
1091in the \textit{namzdf} namelist.
1092
1093This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the
1094bottom boundary condition is implemented implicitly.
1095
1096\begin{equation} \label{Eq_dynzdf_bfr}
1097\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
1098    = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
1099\end{equation}
1100
1101where $mbk$ is the layer number of the bottom wet layer. superscript $n+1$ means the velocity used in the
1102friction formula is to be calculated, so, it is implicit.
1103
1104If split-explicit time splitting is used, care must be taken to avoid the double counting of
1105the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic
1106pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove
1107the bottom friction induced by these two terms which has been included in the 3-D momentum trend
1108and update it with the latest value. On the other hand, the bottom friction contributed by the
1109other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations
1110and should not be added in the 2-D barotropic mode.
1111
1112The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the
1113following:
1114
1115\begin{equation} \label{Eq_dynspg_ts_bfr1}
1116\frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
1117\left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
1118\end{equation}
1119\begin{equation} \label{Eq_dynspg_ts_bfr2}
1120\frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
1121\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
11222\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
1123\end{equation}
1124
1125where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping
1126is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
1127 $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops
1128while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom
1129layer horizontal velocity.
1130
1131
1132
1133
1134% -------------------------------------------------------------------------------------------------------------
1135%       Bottom Friction with split-explicit time splitting
1136% -------------------------------------------------------------------------------------------------------------
1137\subsection[Bottom Friction with split-explicit time splitting]{Bottom Friction with split-explicit time splitting (\protect\np{ln\_bfrimp})}
1138\label{ZDF_bfr_ts}
1139
1140When calculating the momentum trend due to bottom friction in \mdl{dynbfr}, the
1141bottom velocity at the before time step is used. This velocity includes both the
1142baroclinic and barotropic components which is appropriate when using either the
1143explicit or filtered surface pressure gradient algorithms (\key{dynspg\_exp} or
1144\key{dynspg\_flt}). Extra attention is required, however, when using
1145split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface
1146equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three
1147dimensional prognostic variables are solved with the longer time step
1148of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom
1149friction appropriate to this method is that given by the selected parameterisation
1150($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities
1151at each barotropic timestep.
1152
1153In the case of non-linear bottom friction, we have elected to partially linearise
1154the problem by keeping the coefficients fixed throughout the barotropic
1155time-stepping to those computed in \mdl{zdfbfr} using the now timestep.
1156This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to:
1157
1158\begin{enumerate}
1159\item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before
1160barotropic velocity to the bottom friction component of the vertically
1161integrated momentum trend. Note the same stability check that is carried out
1162on the bottom friction coefficient in \mdl{dynbfr} has to be applied here to
1163ensure that the trend removed matches that which was added in \mdl{dynbfr}.
1164\item At each barotropic step, compute the contribution of the current barotropic
1165velocity to the trend due to bottom friction. Add this contribution to the
1166vertically integrated momentum trend. This contribution is handled implicitly which
1167eliminates the need to impose a stability criteria on the values of the bottom friction
1168coefficient within the barotropic loop.
1169\end{enumerate}
1170
1171Note that the use of an implicit formulation within the barotropic loop
1172for the bottom friction trend means that any limiting of the bottom friction coefficient
1173in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time
1174splitting. This is because the major contribution to bottom friction is likely to come from
1175the barotropic component which uses the unrestricted value of the coefficient. However, if the
1176limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas
1177applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} )
1178which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}.
1179
1180Otherwise, the implicit formulation takes the form:
1181\begin{equation} \label{Eq_zdfbfr_implicitts}
1182 \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] 
1183\end{equation}
1184where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height),
1185$c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and $RHS$ represents
1186all the components to the vertically integrated momentum trend except for that due to bottom friction.
1187
1188
1189
1190
1191% ================================================================
1192% Tidal Mixing
1193% ================================================================
1194\section{Tidal Mixing (\protect\key{zdftmx})}
1195\label{ZDF_tmx}
1196
1197%--------------------------------------------namzdf_tmx--------------------------------------------------
1198%\forfile{../namelists/namzdf_tmx}
1199%--------------------------------------------------------------------------------------------------------------
1200
1201
1202% -------------------------------------------------------------------------------------------------------------
1203%        Bottom intensified tidal mixing
1204% -------------------------------------------------------------------------------------------------------------
1205\subsection{Bottom intensified tidal mixing}
1206\label{ZDF_tmx_bottom}
1207
1208Options are defined through the  \ngn{namzdf\_tmx} namelist variables.
1209The parameterization of tidal mixing follows the general formulation for
1210the vertical eddy diffusivity proposed by \citet{St_Laurent_al_GRL02} and
1211first introduced in an OGCM by \citep{Simmons_al_OM04}.
1212In this formulation an additional vertical diffusivity resulting from internal tide breaking,
1213$A^{vT}_{tides}$ is expressed as a function of $E(x,y)$, the energy transfer from barotropic
1214tides to baroclinic tides :
1215\begin{equation} \label{Eq_Ktides}
1216A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }
1217\end{equation}
1218where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency
1219(see \S\ref{TRA_bn2}), $\rho$ the density, $q$ the tidal dissipation efficiency,
1220and $F(z)$ the vertical structure function.
1221
1222The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter)
1223and is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980).
1224The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter)
1225represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally,
1226with the remaining $1-q$ radiating away as low mode internal waves and
1227contributing to the background internal wave field. A value of $q=1/3$ is typically used 
1228\citet{St_Laurent_al_GRL02}.
1229The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical.
1230It is implemented as a simple exponential decaying upward away from the bottom,
1231with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter, with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04},
1232\begin{equation} \label{Eq_Fz}
1233F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
1234\end{equation}
1235and is normalized so that vertical integral over the water column is unity.
1236
1237The associated vertical viscosity is calculated from the vertical
1238diffusivity assuming a Prandtl number of 1, $i.e.$ $A^{vm}_{tides}=A^{vT}_{tides}$.
1239In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity
1240is capped at $300\,cm^2/s$ and impose a lower limit on $N^2$ of \np{rn\_n2min} 
1241usually set to $10^{-8} s^{-2}$. These bounds are usually rarely encountered.
1242
1243The internal wave energy map, $E(x, y)$ in \eqref{Eq_Ktides}, is derived
1244from a barotropic model of the tides utilizing a parameterization of the
1245conversion of barotropic tidal energy into internal waves.
1246The essential goal of the parameterization is to represent the momentum
1247exchange between the barotropic tides and the unrepresented internal waves
1248induced by the tidal flow over rough topography in a stratified ocean.
1249In the current version of \NEMO, the map is built from the output of
1250the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.
1251This model provides the dissipation associated with internal wave energy for the M2 and K1
1252tides component (Fig.~\ref{Fig_ZDF_M2_K1_tmx}). The S2 dissipation is simply approximated
1253as being $1/4$ of the M2 one. The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$.
1254Its global mean value is $1.1$ TW, in agreement with independent estimates
1255\citep{Egbert_Ray_Nat00, Egbert_Ray_JGR01}.
1256
1257%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1258\begin{figure}[!t]   \begin{center}
1259\includegraphics[width=0.90\textwidth]{Fig_ZDF_M2_K1_tmx}
1260\caption{  \protect\label{Fig_ZDF_M2_K1_tmx} 
1261(a) M2 and (b) K1 internal wave drag energy from \citet{Carrere_Lyard_GRL03} ($W/m^2$). }
1262\end{center}   \end{figure}
1263%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1264 
1265% -------------------------------------------------------------------------------------------------------------
1266%        Indonesian area specific treatment
1267% -------------------------------------------------------------------------------------------------------------
1268\subsection{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}
1269\label{ZDF_tmx_itf}
1270
1271When the Indonesian Through Flow (ITF) area is included in the model domain,
1272a specific treatment of tidal induced mixing in this area can be used.
1273It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide
1274an input NetCDF file, \ifile{mask\_itf}, which contains a mask array defining the ITF area
1275where the specific treatment is applied.
1276
1277When \np{ln\_tmx\_itf}=true, the two key parameters $q$ and $F(z)$ are adjusted following
1278the parameterisation developed by \citet{Koch-Larrouy_al_GRL07}:
1279
1280First, the Indonesian archipelago is a complex geographic region
1281with a series of large, deep, semi-enclosed basins connected via
1282numerous narrow straits. Once generated, internal tides remain
1283confined within this semi-enclosed area and hardly radiate away.
1284Therefore all the internal tides energy is consumed within this area.
1285So it is assumed that $q = 1$, $i.e.$ all the energy generated is available for mixing.
1286Note that for test purposed, the ITF tidal dissipation efficiency is a
1287namelist parameter (\np{rn\_tfe\_itf}). A value of $1$ or close to is
1288this recommended for this parameter.
1289
1290Second, the vertical structure function, $F(z)$, is no more associated
1291with a bottom intensification of the mixing, but with a maximum of
1292energy available within the thermocline. \citet{Koch-Larrouy_al_GRL07} 
1293have suggested that the vertical distribution of the energy dissipation
1294proportional to $N^2$ below the core of the thermocline and to $N$ above.
1295The resulting $F(z)$ is:
1296\begin{equation} \label{Eq_Fz_itf}
1297F(i,j,k) \sim     \left\{ \begin{aligned}
1298\frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\
1299\frac{q\,\Gamma E(i,j) } {\rho     \, \int N^2 dz}    \qquad \text{when $\partial_z N > 0$}
1300                      \end{aligned} \right.
1301\end{equation}
1302
1303Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$,
1304which agrees with the independent estimates inferred from observations.
1305Introduced in a regional OGCM, the parameterization improves the water mass
1306characteristics in the different Indonesian seas, suggesting that the horizontal
1307and vertical distributions of the mixing are adequately prescribed
1308\citep{Koch-Larrouy_al_GRL07, Koch-Larrouy_al_OD08a, Koch-Larrouy_al_OD08b}.
1309Note also that such a parameterisation has a significant impact on the behaviour
1310of global coupled GCMs \citep{Koch-Larrouy_al_CD10}.
1311
1312
1313% ================================================================
1314% Internal wave-driven mixing
1315% ================================================================
1316\section{Internal wave-driven mixing (\protect\key{zdftmx\_new})}
1317\label{ZDF_tmx_new}
1318
1319%--------------------------------------------namzdf_tmx_new------------------------------------------
1320%\forfile{../namelists/namzdf_tmx_new}
1321%--------------------------------------------------------------------------------------------------------------
1322
1323The parameterization of mixing induced by breaking internal waves is a generalization
1324of the approach originally proposed by \citet{St_Laurent_al_GRL02}.
1325A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
1326and the resulting diffusivity is obtained as
1327\begin{equation} \label{Eq_Kwave}
1328A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
1329\end{equation}
1330where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution
1331of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false,
1332the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}.
1333In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter
1334$Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater,
1335following the model of \cite{Bouffard_Boegman_DAO2013} 
1336and the implementation of \cite{de_lavergne_JPO2016_efficiency}.
1337Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant.
1338
1339In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
1340as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice).
1341This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014},
1342is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
1343
1344The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed
1345from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$,
1346$E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures
1347(de Lavergne et al., in prep):
1348\begin{align*}
1349F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
1350F_{pyc}(i,j,k) &\propto N^{n\_p}\\
1351F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
1352\end{align*} 
1353In the above formula, $h_{ab}$ denotes the height above bottom,
1354$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
1355\begin{equation*}
1356h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
1357\end{equation*}
1358The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation.
1359It can take values of 1 (recommended) or 2.
1360Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of
1361the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.
1362$h_{cri}$ is related to the large-scale topography of the ocean (etopo2)
1363and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
1364the abyssal hill topography \citep{Goff_JGR2010} and the latitude.
1365
1366% ================================================================
1367
1368
1369
1370\end{document}
Note: See TracBrowser for help on using the repository browser.