source: trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex @ 6497

Last change on this file since 6497 was 6497, checked in by gm, 5 years ago

#1720 - trunk: add Casimir tidal parameterization

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