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

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

Reinstate svn keywords

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