source: branches/UKMO/v3_6_extra_CMIP6_diagnostics/DOC/TexFiles/Chapters/Chap_ZDF.tex @ 7492

Last change on this file since 7492 was 7492, checked in by timgraham, 5 years ago

Final merge of head of v3.6 STABLE branch into this branch

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