Context Navigation

← Previous Changeset
Next Changeset →

Changeset 11213

Timestamp:

2019-07-04T11:29:44+02:00 (5 years ago)

Author:

jchanut

Message:

#2216: Add Romain's update of ZDF chapter

Location:

NEMO/trunk/doc/latex/NEMO

Files:

: 2 edited

main/bibliography.bib (modified) (2 diffs)
subfiles/chap_ZDF.tex (modified) (45 diffs)

Legend:

: Unmodified
: Added
: Removed

NEMO/trunk/doc/latex/NEMO/main/bibliography.bib

-                      r11170
+                      r11213
+}
+@article{         reffray.guillaume.ea_GMD15,
+  title         = "Modelling turbulent vertical mixing sensitivity using
+                   a 1-D version of NEMO",
+  pages         = "69--86",
+  journal       = "Geoscientific Model Development",
+  volume        = "8",
+  number        = "1",
+  author        = "G. Reffray. and R. Bourdalle-Badie and C. Calone",
+  year          = "2015",
+  month         = "jan",
+  publisher     = "Copernicus GmbH",
+  issn          = "1991-9603",
+  doi           = "10.5194/gmd-8-69-2015"
+}
 @book{            richtmyer.morton_bk67,
   pages         = "405",
 …
+}
+@article{         shchepetkin_OM15,
+  title         = "An Adaptive, Courant-number-dependent implicit
+       scheme for vertical advection in oceanic modeling",
+  pages         = "38--69",
+  journal       = "Ocean Modelling",
+  volume        = "91",
+  author        = "A. F. Shchepetkin",
+  year          = "2015",
+  month         = "july",
+  publisher     = "Elsevier BV",
+  issn          = "1463-5003",
+  doi           = "10.1016/j.ocemod.2015.03.006"
+}
 @article{         siddorn.furner_OM13,
   title         = "An analytical stretching function that combines the best

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

-                      r11179
+                      r11213
 At the surface they are prescribed from the surface forcing (see \autoref{chap:SBC}),
 while at the bottom they are set to zero for heat and salt,
 unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \key{trabbl} defined,
+unless a geothermal flux forcing is prescribed as a bottom boundary condition (\ie \np{ln\_trabbc} defined,
 see \autoref{subsec:TRA_bbc}), and specified through a bottom friction parameterisation for momentum
 (see \autoref{sec:ZDF_bfr}).
+(see \autoref{sec:ZDF_drg}).
 In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and
 …
 respectively (see \autoref{sec:TRA_zdf} and \autoref{sec:DYN_zdf}).
 These coefficients can be assumed to be either constant, or a function of the local Richardson number,
 or computed from a turbulent closure model (either TKE or GLS formulation).
 The computation of these coefficients is initialized in the \mdl{zdfini} module and performed in
 the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} modules.
+or computed from a turbulent closure model (either TKE or GLS or OSMOSIS formulation).
+The computation of these coefficients is initialized in the \mdl{zdfphy} module and performed in
+the \mdl{zdfric}, \mdl{zdftke} or \mdl{zdfgls} or \mdl{zdfosm} modules.
 The trends due to the vertical momentum and tracer diffusion, including the surface forcing,
 are computed and added to the general trend in the \mdl{dynzdf} and \mdl{trazdf} modules, respectively.
+These trends can be computed using either a forward time stepping scheme
+(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme
+(\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients,
+and thus of the formulation used (see \autoref{chap:STP}).
+% -------------------------------------------------------------------------------------------------------------
+%        Constant
+% -------------------------------------------------------------------------------------------------------------
+\subsection[Constant (\texttt{\textbf{key\_zdfcst}})]
+{Constant (\protect\key{zdfcst})}
+\label{subsec:ZDF_cst}
+%--------------------------------------------namzdf---------------------------------------------------------
+%These trends can be computed using either a forward time stepping scheme
+%(namelist parameter \np{ln\_zdfexp}\forcode{ = .true.}) or a backward time stepping scheme
+%(\np{ln\_zdfexp}\forcode{ = .false.}) depending on the magnitude of the mixing coefficients,
+%and thus of the formulation used (see \autoref{chap:STP}).
+%--------------------------------------------namzdf--------------------------------------------------------
 \nlst{namzdf}
 %--------------------------------------------------------------------------------------------------------------
+% -------------------------------------------------------------------------------------------------------------
+%        Constant
+% -------------------------------------------------------------------------------------------------------------
+\subsection[Constant (\forcode{ln_zdfcst = .true.})]
+{Constant (\protect\np{ln\_zdfcst}\forcode{ = .true.})}
+\label{subsec:ZDF_cst}
 Options are defined through the \ngn{namzdf} namelist variables.
 When \key{zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to
+When \np{ln\_zdfcst} is defined, the momentum and tracer vertical eddy coefficients are set to
 constant values over the whole ocean.
 This is the crudest way to define the vertical ocean physics.
 It is recommended that this option is only used in process studies, not in basin scale simulations.
+It is recommended to use this option only in process studies, not in basin scale simulations.
 Typical values used in this case are:
 \begin{align*}
 …
 %        Richardson Number Dependent
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Richardson number dependent (\texttt{\textbf{key\_zdfric}})]
 {Richardson number dependent (\protect\key{zdfric})}
+\subsection[Richardson number dependent (\forcode{ln_zdfric = .true.})]
+{Richardson number dependent (\protect\np{ln\_zdfric}\forcode{ = .true.})}
 \label{subsec:ZDF_ric}
 …
 %--------------------------------------------------------------------------------------------------------------
 When \key{zdfric} is defined, a local Richardson number dependent formulation for the vertical momentum and
+When \np{ln\_zdfric}\forcode{ = .true.}, a local Richardson number dependent formulation for the vertical momentum and
 tracer eddy coefficients is set through the \ngn{namzdf\_ric} namelist variables.
 The vertical mixing coefficients are diagnosed from the large scale variables computed by the model.
 …
 %        TKE Turbulent Closure Scheme
 % -------------------------------------------------------------------------------------------------------------
 \subsection[TKE turbulent closure scheme (\texttt{\textbf{key\_zdftke}})]
 {TKE turbulent closure scheme (\protect\key{zdftke})}
+\subsection[TKE turbulent closure scheme (\forcode{ln_zdftke = .true.})]
+{TKE turbulent closure scheme (\protect\np{ln\_zdftke}\forcode{ = .true.})}
 \label{subsec:ZDF_tke}
 %--------------------------------------------namzdf_tke--------------------------------------------------
 …
   \end{cases}
 \end{align*}
-Options are defined through the  \ngn{namzdfy\_tke} namelist variables.
 The choice of $P_{rt}$ is controlled by the \np{nn\_pdl} namelist variable.
 …
 too weak vertical diffusion.
 They must be specified at least larger than the molecular values, and are set through \np{rn\_avm0} and
 \np{rn\_avt0} (namzdf namelist, see \autoref{subsec:ZDF_cst}).
+\np{rn\_avt0} (\ngn{namzdf} namelist, see \autoref{subsec:ZDF_cst}).
 \subsubsection{Turbulent length scale}
 …
 it makes no sense for locally unstable stratification and the computation no longer uses all
 the information contained in the vertical density profile.
 To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2..3} cases,
+To overcome these drawbacks, \citet{madec.delecluse.ea_NPM98} introduces the \np{nn\_mxl}\forcode{ = 2, 3} cases,
 which add an extra assumption concerning the vertical gradient of the computed length scale.
 So, the length scales are first evaluated as in \autoref{eq:tke_mxl0_1} and then bounded such that:
 …
 The \citet{mellor.blumberg_JPO04} modifications acts on surface length scale and TKE values and
 air-sea drag coefficient.
 The latter concerns the bulk formulea and is not discussed here.
+The latter concerns the bulk formulae and is not discussed here.
 Following \citet{craig.banner_JPO94}, the boundary condition on surface TKE value is :
 …
 with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}\forcode{ = 67.83} corresponds
 to $\alpha_{CB} = 100$.
 Further setting  \np{ln\_mxl0} to true applies \autoref{eq:ZDF_Lsbc} as surface boundary condition on length scale,
+Further setting  \np{ln\_mxl0=.true.},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,
 with $\beta$ hard coded to the Stacey's value.
 Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on
+Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the
 surface $\bar{e}$ value.
 …
 \citep{axell_JGR02} for a $k-\epsilon$ turbulent closure.
 The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in
 an extra source terms of TKE, $P_{LC}$.
+an extra source term of TKE, $P_{LC}$.
 The presence of $P_{LC}$ in \autoref{eq:zdftke_e}, the TKE equation, is controlled by setting \np{ln\_lc} to
 \forcode{.true.} in the namtke namelist.
+\forcode{.true.} in the \ngn{namzdf\_tke} namelist.
 By making an analogy with the characteristic convective velocity scale (\eg, \citet{dalessio.abdella.ea_JPO98}),
 …
 The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations:
 $H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by
+$H_{LC}$ is the depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by
 converting its kinetic energy to potential energy, according to
 \[
 …
 \end{equation}
 where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that
 penetrate in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of
+penetrates in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of
 the penetration, and $f_i$ is the ice concentration
 (no penetration if $f_i=1$, that is if the ocean is entirely covered by sea-ice).
+(no penetration if $f_i=1$, \ie if the ocean is entirely covered by sea-ice).
 The value of $f_r$, usually a few percents, is specified through \np{rn\_efr} namelist parameter.
 The vertical mixing length scale, $h_\tau$, can be set as a 10~m uniform value (\np{nn\_etau}\forcode{ = 0}) or
 …
 (\np{nn\_etau}\forcode{ = 1}).
 Note that two other option existe, \np{nn\_etau}\forcode{ = 2..3}.
+Note that two other option exist, \np{nn\_etau}\forcode{ = 2, 3}.
 They correspond to applying \autoref{eq:ZDF_Ehtau} only at the base of the mixed layer,
 or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrate the ocean.
+or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrates the ocean.
 Those two options are obsolescent features introduced for test purposes.
 They will be removed in the next release.
+% This should be explain better below what this rn_eice parameter is meant for:
+In presence of Sea Ice, the value of this mixing can be modulated by the \np{rn\_eice} namelist parameter.
+This parameter varies from \forcode{0} for no effect to \forcode{4} to suppress the TKE input into the ocean when Sea Ice concentration
+is greater than 25\%.
 % from Burchard et al OM 2008 :
 …
 % -------------------------------------------------------------------------------------------------------------
-%        TKE discretization considerations
-% -------------------------------------------------------------------------------------------------------------
-\subsection[TKE discretization considerations (\texttt{\textbf{key\_zdftke}})]
-{TKE discretization considerations (\protect\key{zdftke})}
-\label{subsec:ZDF_tke_ene}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]
-  \begin{center}
-    \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme}
-    \caption{
-      \protect\label{fig:TKE_time_scheme}
-      Illustration of the TKE time integration and its links to the momentum and tracer time integration.
+    }
-  \end{center}
-\end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-The production of turbulence by vertical shear (the first term of the right hand side of
-\autoref{eq:zdftke_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion
-(first line in \autoref{eq:PE_zdf}).
-To do so a special care have to be taken for both the time and space discretization of
-the TKE equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}.
-Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how
-the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with
-the one-level forward time stepping of TKE equation.
-With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to
-the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and
-summing the result vertically:
-\begin{equation}
-  \label{eq:energ1}
-  \begin{split}
-    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
-    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}
-    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
-  \end{split}
-\end{equation}
-Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$,
-known at time $t$ (\autoref{fig:TKE_time_scheme}), as it is required when using the TKE scheme
-(see \autoref{sec:STP_forward_imp}).
-The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy transfer at
-the surface (atmospheric forcing) and at the bottom (friction effect).
-The second term is always negative.
-It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$.
-\autoref{eq:energ1} implies that, to be energetically consistent,
-the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
-${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$
-(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
-A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
-(second term of the right hand side of \autoref{eq:zdftke_e}).
-This term must balance the input of potential energy resulting from vertical mixing.
-The rate of change of potential energy (in 1D for the demonstration) due vertical mixing is obtained by
-multiplying vertical density diffusion tendency by $g\,z$ and and summing the result vertically:
-\begin{equation}
-  \label{eq:energ2}
-  \begin{split}
-    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
-    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}
-    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
-    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
-    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
-  \end{split}
-\end{equation}
-where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
-The first term of the right hand side of \autoref{eq:energ2} is always zero because
-there is no diffusive flux through the ocean surface and bottom).
-The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
-Therefore \autoref{eq:energ1} implies that, to be energetically consistent,
-the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e}, the TKE equation.
-Let us now address the space discretization issue.
-The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in
-the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:cell}).
-A space averaging is thus required to obtain the shear TKE production term.
-By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by
-the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
-Furthermore, the possible time variation of $e_3$ (\key{vvl} case) have to be taken into account.
-The above energetic considerations leads to the following final discrete form for the TKE equation:
-\begin{equation}
-  \label{eq:zdftke_ene}
-  \begin{split}
-    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv
-    \Biggl\{ \Biggr.
-    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} }
-        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
-    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} }
-        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j}
-    \Biggr. \Biggr\}   \\
+    %
-    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
+    %
-    +&\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]   \\
+    %
-    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
-  \end{split}
-\end{equation}
-where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation)
-are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}).
-Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible.
-The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as
-they all appear in the right hand side of \autoref{eq:zdftke_ene}.
-For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
-% -------------------------------------------------------------------------------------------------------------
 %        GLS Generic Length Scale Scheme
 % -------------------------------------------------------------------------------------------------------------
 \subsection[GLS: Generic Length Scale (\texttt{\textbf{key\_zdfgls}})]
 {GLS: Generic Length Scale (\protect\key{zdfgls})}
+\subsection[GLS: Generic Length Scale (\forcode{ln_zdfgls = .true.})]
+{GLS: Generic Length Scale (\protect\np{ln\_zdfgls}\forcode{ = .true.})}
 \label{subsec:ZDF_gls}
 …
 The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of
 the choice of the turbulence model.
 Four different turbulent models are pre-defined (Tab.\autoref{tab:GLS}).
+Four different turbulent models are pre-defined (\autoref{tab:GLS}).
 They are made available through the \np{nn\_clo} namelist parameter.
 …
       \protect\label{tab:GLS}
       Set of predefined GLS parameters, or equivalently predefined turbulence models available with
       \protect\key{zdfgls} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}.
+      \protect\np{ln\_zdfgls}\forcode{ = .true.} and controlled by the \protect\np{nn\_clos} namelist variable in \protect\ngn{namzdf\_gls}.
+    }
   \end{center}
 …
 In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of
 the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries
+the mixing length towards $\kappa z_b$ ($\kappa$ is the Von Karman constant and $z_b$ the rugosity length scale) value near physical boundaries
 (logarithmic boundary layer law).
 $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by \citet{galperin.kantha.ea_JAS88},
 or by \citet{kantha.clayson_JGR94} or one of the two functions suggested by \citet{canuto.howard.ea_JPO01}
 (\np{nn\_stab\_func}\forcode{ = 0..3}, resp.).
 The value of $C_{0\mu}$ depends of the choice of the stability function.
+(\np{nn\_stab\_func}\forcode{ = 0, 3}, resp.).
+The value of $C_{0\mu}$ depends on the choice of the stability function.
 The surface and bottom boundary condition on both $\bar{e}$ and $\psi$ can be calculated thanks to Dirichlet or
 …
 The time and space discretization of the GLS equations follows the same energetic consideration as for
 the TKE case described in \autoref{subsec:ZDF_tke_ene} \citep{burchard_OM02}.
+Examples of performance of the 4 turbulent closure scheme can be found in \citet{warner.sherwood.ea_OM05}.
+Evaluation of the 4 GLS turbulent closure schemes can be found in \citet{warner.sherwood.ea_OM05} in ROMS model and
+ in \citet{reffray.guillaume.ea_GMD15} for the \NEMO model.
 % -------------------------------------------------------------------------------------------------------------
 %        OSM OSMOSIS BL Scheme
 % -------------------------------------------------------------------------------------------------------------
 \subsection[OSM: OSMosis boundary layer scheme (\texttt{\textbf{key\_zdfosm}})]
 {OSM: OSMosis boundary layer scheme (\protect\key{zdfosm})}
+\subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm = .true.})]
+{OSM: OSMosis boundary layer scheme (\protect\np{ln\_zdfosm}\forcode{ = .true.})}
 \label{subsec:ZDF_osm}
 %--------------------------------------------namzdf_osm---------------------------------------------------------
 …
 The OSMOSIS turbulent closure scheme is based on......   TBC
+% -------------------------------------------------------------------------------------------------------------
+%        TKE and GLS discretization considerations
+% -------------------------------------------------------------------------------------------------------------
+\subsection[ Discrete energy conservation for TKE and GLS schemes]
+{Discrete energy conservation for TKE and GLS schemes}
+\label{subsec:ZDF_tke_ene}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_ZDF_TKE_time_scheme}
+    \caption{
+      \protect\label{fig:TKE_time_scheme}
+      Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and its links to the momentum and tracer time integration.
+    }
+  \end{center}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+The production of turbulence by vertical shear (the first term of the right hand side of
+\autoref{eq:zdftke_e}) and  \autoref{eq:zdfgls_e}) should balance the loss of kinetic energy associated with the vertical momentum diffusion
+(first line in \autoref{eq:PE_zdf}).
+To do so a special care has to be taken for both the time and space discretization of
+the kinetic energy equation \citep{burchard_OM02,marsaleix.auclair.ea_OM08}.
+Let us first address the time stepping issue. \autoref{fig:TKE_time_scheme} shows how
+the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with
+the one-level forward time stepping of the equation for $\bar{e}$.
+With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to
+the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and
+summing the result vertically:
+\begin{equation}
+  \label{eq:energ1}
+  \begin{split}
+    \int_{-H}^{\eta}  u^t \,\partial_z &\left( {K_m}^t \,(\partial_z u)^{t+\rdt}  \right) \,dz   \\
+    &= \Bigl[  u^t \,{K_m}^t \,(\partial_z u)^{t+\rdt} \Bigr]_{-H}^{\eta}
+    - \int_{-H}^{\eta}{ {K_m}^t \,\partial_z{u^t} \,\partial_z u^{t+\rdt} \,dz }
+  \end{split}
+\end{equation}
+Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$,
+known at time $t$ (\autoref{fig:TKE_time_scheme}), as it is required when using the TKE scheme
+(see \autoref{sec:STP_forward_imp}).
+The first term of the right hand side of \autoref{eq:energ1} represents the kinetic energy transfer at
+the surface (atmospheric forcing) and at the bottom (friction effect).
+The second term is always negative.
+It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar{e}$.
+\autoref{eq:energ1} implies that, to be energetically consistent,
+the production rate of $\bar{e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as
+${K_m}^{t-\rdt}\,(\partial_z u)^{t-\rdt} \,(\partial_z u)^t$
+(and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).
+A similar consideration applies on the destruction rate of $\bar{e}$ due to stratification
+(second term of the right hand side of \autoref{eq:zdftke_e} and \autoref{eq:zdfgls_e}).
+This term must balance the input of potential energy resulting from vertical mixing.
+The rate of change of potential energy (in 1D for the demonstration) due to vertical mixing is obtained by
+multiplying the vertical density diffusion tendency by $g\,z$ and and summing the result vertically:
+\begin{equation}
+  \label{eq:energ2}
+  \begin{split}
+    \int_{-H}^{\eta} g\,z\,\partial_z &\left( {K_\rho}^t \,(\partial_k \rho)^{t+\rdt}   \right) \,dz    \\
+    &= \Bigl[  g\,z \,{K_\rho}^t \,(\partial_z \rho)^{t+\rdt} \Bigr]_{-H}^{\eta}
+    - \int_{-H}^{\eta}{ g \,{K_\rho}^t \,(\partial_k \rho)^{t+\rdt} } \,dz   \\
+    &= - \Bigl[  z\,{K_\rho}^t \,(N^2)^{t+\rdt} \Bigr]_{-H}^{\eta}
+    + \int_{-H}^{\eta}{  \rho^{t+\rdt} \, {K_\rho}^t \,(N^2)^{t+\rdt} \,dz  }
+  \end{split}
+\end{equation}
+where we use $N^2 = -g \,\partial_k \rho / (e_3 \rho)$.
+The first term of the right hand side of \autoref{eq:energ2} is always zero because
+there is no diffusive flux through the ocean surface and bottom).
+The second term is minus the destruction rate of  $\bar{e}$ due to stratification.
+Therefore \autoref{eq:energ1} implies that, to be energetically consistent,
+the product ${K_\rho}^{t-\rdt}\,(N^2)^t$ should be used in \autoref{eq:zdftke_e} and  \autoref{eq:zdfgls_e}.
+Let us now address the space discretization issue.
+The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in
+the centre of the side faces of a $t$-box in staggered C-grid (\autoref{fig:cell}).
+A space averaging is thus required to obtain the shear TKE production term.
+By redoing the \autoref{eq:energ1} in the 3D case, it can be shown that the product of eddy coefficient by
+the shear at $t$ and $t-\rdt$ must be performed prior to the averaging.
+Furthermore, the time variation of $e_3$ has be taken into account.
+The above energetic considerations leads to the following final discrete form for the TKE equation:
+\begin{equation}
+  \label{eq:zdftke_ene}
+  \begin{split}
+    \frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt}  \equiv
+    \Biggl\{ \Biggr.
+    &\overline{ \left( \left(\overline{K_m}^{\,i+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[u^{t+\rdt}]}{{e_3u}^{t+\rdt} }
+        \ \frac{\delta_{k+1/2}[u^ t         ]}{{e_3u}^ t          }  \right) }^{\,i} \\
+    +&\overline{  \left( \left(\overline{K_m}^{\,j+1/2}\right)^{t-\rdt} \,\frac{\delta_{k+1/2}[v^{t+\rdt}]}{{e_3v}^{t+\rdt} }
+        \ \frac{\delta_{k+1/2}[v^ t         ]}{{e_3v}^ t          }  \right) }^{\,j}
+    \Biggr. \Biggr\}   \\
+    %
+    - &{K_\rho}^{t-\rdt}\,{(N^2)^t}    \\
+    %
+    +&\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]   \\
+    %
+    - &c_\epsilon \; \left( \frac{\sqrt{\bar {e}}}{l_\epsilon}\right)^{t-\rdt}\,(\bar {e})^{t+\rdt}
+  \end{split}
+\end{equation}
+where the last two terms in \autoref{eq:zdftke_ene} (vertical diffusion and Kolmogorov dissipation)
+are time stepped using a backward scheme (see\autoref{sec:STP_forward_imp}).
+Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible.
+%The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as
+%they all appear in the right hand side of \autoref{eq:zdftke_ene}.
+%For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.
 % ================================================================
 % Convection
 …
 \section{Convection}
 \label{sec:ZDF_conv}
-%--------------------------------------------namzdf--------------------------------------------------------
-\nlst{namzdf}
-%--------------------------------------------------------------------------------------------------------------
 Static instabilities (\ie light potential densities under heavy ones) may occur at particular ocean grid points.
 …
 {Non-penetrative convective adjustment (\protect\np{ln\_tranpc}\forcode{ = .true.})}
 \label{subsec:ZDF_npc}
-%--------------------------------------------namzdf--------------------------------------------------------
-\nlst{namzdf}
-%--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 Two main differences have been introduced compared to the original algorithm:
 $(i)$ the stability is now checked using the Brunt-V\"{a}is\"{a}l\"{a} frequency
 (not the the difference in potential density);
+(not the difference in potential density);
 $(ii)$ when two levels are found unstable, their thermal and haline expansion coefficients are vertically mixed in
 the same way their temperature and salinity has been mixed.
 …
 \label{subsec:ZDF_evd}
-%--------------------------------------------namzdf--------------------------------------------------------
-\nlst{namzdf}
-%--------------------------------------------------------------------------------------------------------------
 Options are defined through the  \ngn{namzdf} namelist variables.
 The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}\forcode{ = .true.}.
 In this case, the vertical eddy mixing coefficients are assigned very large values
 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable
+In this case, the vertical eddy mixing coefficients are assigned very large values
+in regions where the stratification is unstable
 (\ie when $N^2$ the Brunt-Vais\"{a}l\"{a} frequency is negative) \citep{lazar_phd97, lazar.madec.ea_JPO99}.
 This is done either on tracers only (\np{nn\_evdm}\forcode{ = 0}) or
 …
 the convective adjustment algorithm presented above when mixing both tracers and
 momentum in the case of static instabilities.
-It requires the use of an implicit time stepping on vertical diffusion terms
-(\ie np{ln\_zdfexp}\forcode{ = .false.}).
 Note that the stability test is performed on both \textit{before} and \textit{now} values of $N^2$.
 …
 %       Turbulent Closure Scheme
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Turbulent closure scheme (\texttt{\textbf{key\_zdf}}\texttt{\textbf{\{tke,gls,osm\}}})]
 {Turbulent Closure Scheme (\protect\key{zdftke}, \protect\key{zdfgls} or \protect\key{zdfosm})}
+\subsection[Handling convection with turbulent closure schemes (\forcode{ln_zdf/tke/gls/osm = .true.})]
+{Handling convection with turbulent closure schemes (\protect\np{ln\_zdf/tke/gls/osm}\forcode{ = .true.})}
 \label{subsec:ZDF_tcs}
+The turbulent closure scheme presented in \autoref{subsec:ZDF_tke} and \autoref{subsec:ZDF_gls}
+(\key{zdftke} or \key{zdftke} is defined) in theory solves the problem of statically unstable density profiles.
+The turbulent closure schemes presented in \autoref{subsec:ZDF_tke}, \autoref{subsec:ZDF_gls} and
+\autoref{subsec:ZDF_osm} (\ie \np{ln\_zdftke} or \np{ln\_zdftke} or \np{ln\_zdfosm} defined) deal, in theory,
+with statically unstable density profiles.
 In such a case, the term corresponding to the destruction of turbulent kinetic energy through stratification in
 \autoref{eq:zdftke_e} or \autoref{eq:zdfgls_e} becomes a source term, since $N^2$ is negative.
 It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also the four neighbouring $A_u^{vm} {and}\;A_v^{vm}$
 (up to $1\;m^2s^{-1}$).
+It results in large values of $A_T^{vT}$ and  $A_T^{vT}$, and also of the four neighboring values at
+velocity points $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1}$).
 These large values restore the static stability of the water column in a way similar to that of
 the enhanced vertical diffusion parameterisation (\autoref{subsec:ZDF_evd}).
 …
 It can thus be useful to combine the enhanced vertical diffusion with the turbulent closure scheme,
 \ie setting the \np{ln\_zdfnpc} namelist parameter to true and
 defining the turbulent closure CPP key all together.
 The KPP turbulent closure scheme already includes enhanced vertical diffusion in the case of convection,
 as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
 therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the KPP scheme.
+defining the turbulent closure (\np{ln\_zdftke} or \np{ln\_zdfgls} = \forcode{.true.}) all together.
+The OSMOSIS turbulent closure scheme already includes enhanced vertical diffusion in the case of convection,
+%as governed by the variables $bvsqcon$ and $difcon$ found in \mdl{zdfkpp},
+therefore \np{ln\_zdfevd}\forcode{ = .false.} should be used with the OSMOSIS scheme.
 % gm%  + one word on non local flux with KPP scheme trakpp.F90 module...
 …
 % Double Diffusion Mixing
 % ================================================================
+\section[Double diffusion mixing (\texttt{\textbf{key\_zdfddm}})]
+{Double diffusion mixing (\protect\key{zdfddm})}
+\label{sec:ZDF_ddm}
+\section[Double diffusion mixing (\forcode{ln_zdfddm = .true.})]
+{Double diffusion mixing (\protect\np{ln\_zdfddm}\forcode{ = .true.})}
+\label{subsec:ZDF_ddm}
 %-------------------------------------------namzdf_ddm-------------------------------------------------
 …
 %--------------------------------------------------------------------------------------------------------------
+Options are defined through the  \ngn{namzdf\_ddm} namelist variables.
+This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the namelist parameter
+\np{ln\_zdfddm} in \ngn{namzdf}.
 Double diffusion occurs when relatively warm, salty water overlies cooler, fresher water, or vice versa.
 The former condition leads to salt fingering and the latter to diffusive convection.
 …
 it leads to relatively minor changes in circulation but exerts significant regional influences on
 temperature and salinity.
+This parameterisation has been introduced in \mdl{zdfddm} module and is controlled by the \key{zdfddm} CPP key.
 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients
 …
 % Bottom Friction
 % ================================================================
 \section[Bottom and top friction (\textit{zdfbfr.F90})]
 {Bottom and top friction (\protect\mdl{zdfbfr})}
 \label{sec:ZDF_bfr}
+ \section[Bottom and top friction (\textit{zdfdrg.F90})]
+ {Bottom and top friction (\protect\mdl{zdfdrg})}
+ \label{sec:ZDF_drg}
 %--------------------------------------------nambfr--------------------------------------------------------
+%
+%\nlst{nambfr}
+\nlst{namdrg}
+\nlst{namdrg_top}
+\nlst{namdrg_bot}
 %--------------------------------------------------------------------------------------------------------------
 Options to define the top and bottom friction are defined through the \ngn{nambfr} namelist variables.
+Options to define the top and bottom friction are defined through the \ngn{namdrg} namelist variables.
 The bottom friction represents the friction generated by the bathymetry.
 The top friction represents the friction generated by the ice shelf/ocean interface.
 As the friction processes at the top and bottom are treated in similar way,
 only the bottom friction is described in detail below.
+As the friction processes at the top and the bottom are treated in and identical way,
+the description below considers mostly the bottom friction case, if not stated otherwise.
 …
 a condition on the vertical diffusive flux.
 For the bottom boundary layer, one has:
 \[
   % \label{eq:zdfbfr_flux}
   A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
 \]
+ \[
+   % \label{eq:zdfbfr_flux}
+   A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
+ \]
 where ${\cal F}_h^{\textbf U}$ is represents the downward flux of horizontal momentum outside
 the logarithmic turbulent boundary layer (thickness of the order of 1~m in the ocean).
 …
 To illustrate this, consider the equation for $u$ at $k$, the last ocean level:
 \begin{equation}
   \label{eq:zdfbfr_flux2}
+  \label{eq:zdfdrg_flux2}
   \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}}
 \end{equation}
 …
 In the code, the bottom friction is imposed by adding the trend due to the bottom friction to
 the general momentum trend in \mdl{dynbfr}.
+ the general momentum trend in \mdl{dynzdf}.
 For the time-split surface pressure gradient algorithm, the momentum trend due to
 the barotropic component needs to be handled separately.
 For this purpose it is convenient to compute and store coefficients which can be simply combined with
 bottom velocities and geometric values to provide the momentum trend due to bottom friction.
 These coefficients are computed in \mdl{zdfbfr} and generally take the form $c_b^{\textbf U}$ where:
+ These coefficients are computed in \mdl{zdfdrg} and generally take the form $c_b^{\textbf U}$ where:
 \begin{equation}
   \label{eq:zdfbfr_bdef}
 …
   - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
 \end{equation}
+where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
+where $\textbf{U}_h^b = (u_b\;,\;v_b)$ is the near-bottom, horizontal, ocean velocity.
+Note than from \NEMO 4.0, drag coefficients are only computed at cell centers (\ie at T-points) and refer to as $c_b^T$ in the following. These are then linearly interpolated in space to get $c_b^\textbf{U}$ at velocity points.
 % -------------------------------------------------------------------------------------------------------------
 %       Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Linear bottom friction (\forcode{nn_botfr = [01]})]
 {Linear bottom friction (\protect\np{nn\_botfr}\forcode{ = [01])}}
 \label{subsec:ZDF_bfr_linear}
 The linear bottom friction parameterisation (including the special case of a free-slip condition) assumes that
 the bottom friction is proportional to the interior velocity (\ie the velocity of the last model level):
+ \subsection[Linear top/bottom friction (\forcode{ln_lin = .true.})]
+ {Linear top/bottom friction (\protect\np{ln\_lin}\forcode{ = .true.)}}
+ \label{subsec:ZDF_drg_linear}
+The linear friction parameterisation (including the special case of a free-slip condition) assumes that
+the friction is proportional to the interior velocity (\ie the velocity of the first/last model level):
 \[
   % \label{eq:zdfbfr_linear}
   {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3} \; \frac{\partial \textbf{U}_h}{\partial k} = r \; \textbf{U}_h^b
 \]
 where $r$ is a friction coefficient expressed in ms$^{-1}$.
+where $r$ is a friction coefficient expressed in $m s^{-1}$.
 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean,
 and setting $r = H / \tau$, where $H$ is the ocean depth.
 …
 and assuming an ocean depth $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.
 This is the default value used in \NEMO. It corresponds to a decay time scale of 115~days.
 It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
 For the linear friction case the coefficients defined in the general expression \autoref{eq:zdfbfr_bdef} are:
+It can be changed by specifying \np{rn\_Uc0} (namelist parameter).
+ For the linear friction case the drag coefficient used in the general expression \autoref{eq:zdfbfr_bdef} is:
 \[
   % \label{eq:zdfbfr_linbfr_b}
+  \begin{split}
+    c_b^u &= - r\\
+    c_b^v &= - r\\
+  \end{split}
+\]
+When \np{nn\_botfr}\forcode{ = 1}, the value of $r$ used is \np{rn\_bfri1}.
+Setting \np{nn\_botfr}\forcode{ = 0} is equivalent to setting $r=0$ and
+leads to a free-slip bottom boundary condition.
+These values are assigned in \mdl{zdfbfr}.
+From v3.2 onwards there is support for local enhancement of these values via an externally defined 2D mask array
+(\np{ln\_bfr2d}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file.
+    c_b^T = - r
+\]
+When \forcode{ln_lin = .true.}, the value of $r$ used is \np{rn\_Uc0}*\np{rn\_Cd0}.
+Setting \forcode{ln_OFF = .true.} (and \forcode{ln_lin = .true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
+These values are assigned in \mdl{zdfdrg}.
+Note that there is support for local enhancement of these values via an externally defined 2D mask array
+(\np{ln\_boost}\forcode{ = .true.}) given in the \ifile{bfr\_coef} input NetCDF file.
 The mask values should vary from 0 to 1.
 Locations with a non-zero mask value will have the friction coefficient increased by
 $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
+$mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}.
 % -------------------------------------------------------------------------------------------------------------
 %       Non-Linear Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Non-linear bottom friction (\forcode{nn_botfr = 2})]
 {Non-linear bottom friction (\protect\np{nn\_botfr}\forcode{ = 2})}
 \label{subsec:ZDF_bfr_nonlinear}
 The non-linear bottom friction parameterisation assumes that the bottom friction is quadratic:
 \[
   % \label{eq:zdfbfr_nonlinear}
+ \subsection[Non-linear top/bottom friction (\forcode{ln_no_lin = .true.})]
+ {Non-linear top/bottom friction (\protect\np{ln\_no\_lin}\forcode{ = .true.})}
+ \label{subsec:ZDF_drg_nonlinear}
+The non-linear bottom friction parameterisation assumes that the top/bottom friction is quadratic:
+\[
+  % \label{eq:zdfdrg_nonlinear}
   {\cal F}_h^\textbf{U} = \frac{A^{vm}}{e_3 }\frac{\partial \textbf {U}_h
   }{\partial k}=C_D \;\sqrt {u_b ^2+v_b ^2+e_b } \;\; \textbf {U}_h^b
 \]
 where $C_D$ is a drag coefficient, and $e_b $ a bottom turbulent kinetic energy due to tides,
+where $C_D$ is a drag coefficient, and $e_b $ a top/bottom turbulent kinetic energy due to tides,
 internal waves breaking and other short time scale currents.
 A typical value of the drag coefficient is $C_D = 10^{-3} $.
 …
 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{killworth_JPO92} uses $C_D = 1.4\;10^{-3}$ and
 $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.
+The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} namelist parameters).
+As for the linear case, the bottom friction is imposed in the code by adding the trend due to
+the bottom friction to the general momentum trend in \mdl{dynbfr}.
+For the non-linear friction case the terms computed in \mdl{zdfbfr} are:
+\[
+  % \label{eq:zdfbfr_nonlinbfr}
+  \begin{split}
+    c_b^u &= - \; C_D\;\left[ u^2 + \left(\bar{\bar{v}}^{i+1,j}\right)^2 + e_b \right]^{1/2}\\
+    c_b^v &= - \; C_D\;\left[  \left(\bar{\bar{u}}^{i,j+1}\right)^2 + v^2 + e_b \right]^{1/2}\\
+  \end{split}
+\]
+The coefficients that control the strength of the non-linear bottom friction are initialised as namelist parameters:
+$C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.
+Note for applications which treat tides explicitly a low or even zero value of \np{rn\_bfeb2} is recommended.
+From v3.2 onwards a local enhancement of $C_D$ is possible via an externally defined 2D mask array
+(\np{ln\_bfr2d}\forcode{ = .true.}).
+This works in the same way as for the linear bottom friction case with non-zero masked locations increased by
+$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
+The CME choices have been set as default values (\np{rn\_Cd0} and \np{rn\_ke0} namelist parameters).
+As for the linear case, the friction is imposed in the code by adding the trend due to
+the friction to the general momentum trend in \mdl{dynzdf}.
+For the non-linear friction case the term computed in \mdl{zdfdrg} is:
+\[
+  % \label{eq:zdfdrg_nonlinbfr}
+    c_b^T = - \; C_D\;\left[ \left(\bar{u_b}^{i}\right)^2 + \left(\bar{v_b}^{j}\right)^2 + e_b \right]^{1/2}
+\]
+The coefficients that control the strength of the non-linear friction are initialised as namelist parameters:
+$C_D$= \np{rn\_Cd0}, and $e_b$ =\np{rn\_bfeb2}.
+Note that for applications which consider tides explicitly, a low or even zero value of \np{rn\_bfeb2} is recommended. A local enhancement of $C_D$ is again possible via an externally defined 2D mask array
+(\np{ln\_boost}\forcode{ = .true.}).
+This works in the same way as for the linear friction case with non-zero masked locations increased by
+$mask\_value$ * \np{rn\_boost} * \np{rn\_Cd0}.
 % -------------------------------------------------------------------------------------------------------------
 %       Bottom Friction Log-layer
 % -------------------------------------------------------------------------------------------------------------
+\subsection[Log-layer bottom friction enhancement (\forcode{nn_botfr = 2}, \forcode{ln_loglayer = .true.})]
+{Log-layer bottom friction enhancement (\protect\np{nn\_botfr}\forcode{ = 2}, \protect\np{ln\_loglayer}\forcode{ = .true.})}
+\label{subsec:ZDF_bfr_loglayer}
+In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally enhanced using
+a "law of the wall" scaling.
+If  \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the thickness of
+the last wet layer in each column by:
+\[
+  C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
+\]
+\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness length provided via
+the namelist.
+For stability, the drag coefficient is bounded such that it is kept greater or equal to
+the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional namelist parameter:
+\np{rn\_bfri2\_max}, \ie
+\[
+  rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
+\]
+\noindent Note also that a log-layer enhancement can also be applied to the top boundary friction if
+under ice-shelf cavities are in use (\np{ln\_isfcav}\forcode{ = .true.}).
+In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}.
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction stability
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom friction stability considerations}
+\label{subsec:ZDF_bfr_stability}
+Some care needs to exercised over the choice of parameters to ensure that the implementation of
+bottom friction does not induce numerical instability.
+For the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2} is:
+ \subsection[Log-layer top/bottom friction (\forcode{ln_loglayer = .true.})]
+ {Log-layer top/bottom friction (\protect\np{ln\_loglayer}\forcode{ = .true.})}
+ \label{subsec:ZDF_drg_loglayer}
+In the non-linear friction case, the drag coefficient, $C_D$, can be optionally enhanced using
+a "law of the wall" scaling. This assumes that the model vertical resolution can capture the logarithmic layer which typically occur for layers thinner than 1 m or so.
+If  \np{ln\_loglayer} = .true., $C_D$ is no longer constant but is related to the distance to the wall (or equivalently to the half of the top/bottom layer thickness):
+\[
+  C_D = \left ( {\kappa \over {\mathrm log}\left ( 0.5 \; e_{3b} / rn\_{z0} \right ) } \right )^2
+\]
+\noindent where $\kappa$ is the von-Karman constant and \np{rn\_z0} is a roughness length provided via the namelist.
+The drag coefficient is bounded such that it is kept greater or equal to
+the base \np{rn\_Cd0} value which occurs where layer thicknesses become large and presumably logarithmic layers are not resolved at all. For stability reason, it is also not allowed to exceed the value of an additional namelist parameter:
+\np{rn\_Cdmax}, \ie
+\[
+  rn\_Cd0 \leq C_D \leq rn\_Cdmax
+\]
+\noindent The log-layer enhancement can also be applied to the top boundary friction if
+under ice-shelf cavities are activated (\np{ln\_isfcav}\forcode{ = .true.}).
+%In this case, the relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} and \np{rn\_tfri2\_max}.
+% -------------------------------------------------------------------------------------------------------------
+%       Explicit bottom Friction
+% -------------------------------------------------------------------------------------------------------------
+ \subsection{Explicit top/bottom friction (\forcode{ln_drgimp = .false.})}
+ \label{subsec:ZDF_drg_stability}
+Setting \forcode{ln_drgimp = .false.} means that bottom friction is treated explicitly in time, which has the advantage of simplifying the interaction with the split-explicit free surface (see \autoref{subsec:ZDF_drg_ts}). The latter does indeed require the knowledge of bottom stresses in the course of the barotropic sub-iteration, which becomes less straightforward in the implicit case. In the explicit case, top/bottom stresses can be computed using \textit{before} velocities and inserted in the overall momentum tendency budget. This reads:
+At the top (below an ice shelf cavity):
+\[
+  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t}
+  = c_{t}^{\textbf{U}}\textbf{u}^{n-1}_{t}
+\]
+At the bottom (above the sea floor):
+\[
+  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b}
+  = c_{b}^{\textbf{U}}\textbf{u}^{n-1}_{b}
+\]
+Since this is conditionally stable, some care needs to exercised over the choice of parameters to ensure that the implementation of explicit top/bottom friction does not induce numerical instability.
+For the purposes of stability analysis, an approximation to \autoref{eq:zdfdrg_flux2} is:
 \begin{equation}
   \label{eq:Eqn_bfrstab}
+  \label{eq:Eqn_drgstab}
   \begin{split}
     \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
 …
   \end{split}
 \end{equation}
 \noindent where linear bottom friction and a leapfrog timestep have been assumed.
 To ensure that the bottom friction cannot reverse the direction of flow it is necessary to have:
+\noindent where linear friction and a leapfrog timestep have been assumed.
+To ensure that the friction cannot reverse the direction of flow it is necessary to have:
 \[
   |\Delta u| < \;|u|
 \]
 \noindent which, using \autoref{eq:Eqn_bfrstab}, gives:
+\noindent which, using \autoref{eq:Eqn_drgstab}, gives:
 \[
   r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
 …
 For most applications, with physically sensible parameters these restrictions should not be of concern.
 But caution may be necessary if attempts are made to locally enhance the bottom friction parameters.
 To ensure stability limits are imposed on the bottom friction coefficients both
+To ensure stability limits are imposed on the top/bottom friction coefficients both
 during initialisation and at each time step.
 Checks at initialisation are made in \mdl{zdfbfr} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
+Checks at initialisation are made in \mdl{zdfdrg} (assuming a 1 m.s$^{-1}$ velocity in the non-linear case).
 The number of breaches of the stability criterion are reported as well as
 the minimum and maximum values that have been set.
 The criterion is also checked at each time step, using the actual velocity, in \mdl{dynbfr}.
 Values of the bottom friction coefficient are reduced as necessary to ensure stability;
+The criterion is also checked at each time step, using the actual velocity, in \mdl{dynzdf}.
+Values of the friction coefficient are reduced as necessary to ensure stability;
 these changes are not reported.
 Limits on the bottom friction coefficient are not imposed if the user has elected to
 handle the bottom friction implicitly (see \autoref{subsec:ZDF_bfr_imp}).
+Limits on the top/bottom friction coefficient are not imposed if the user has elected to
+handle the friction implicitly (see \autoref{subsec:ZDF_drg_imp}).
 The number of potential breaches of the explicit stability criterion are still reported for information purposes.
 …
 %       Implicit Bottom Friction
 % -------------------------------------------------------------------------------------------------------------
 \subsection[Implicit bottom friction (\forcode{ln_bfrimp = .true.})]
 {Implicit bottom friction (\protect\np{ln\_bfrimp}\forcode{ = .true.})}
 \label{subsec:ZDF_bfr_imp}
+ \subsection[Implicit top/bottom friction (\forcode{ln_drgimp = .true.})]
+ {Implicit top/bottom friction (\protect\np{ln\_drgimp}\forcode{ = .true.})}
+ \label{subsec:ZDF_drg_imp}
 An optional implicit form of bottom friction has been implemented to improve model stability.
+We recommend this option for shelf sea and coastal ocean applications, especially for split-explicit time splitting.
+This option can be invoked by setting \np{ln\_bfrimp} to \forcode{.true.} in the \textit{nambfr} namelist.
+This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist.
+This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp},
+the bottom boundary condition is implemented implicitly.
+\[
+  % \label{eq:dynzdf_bfr}
+  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
+  = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
+\]
+where $mbk$ is the layer number of the bottom wet layer.
+Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so, it is implicit.
+If split-explicit time splitting is used, care must be taken to avoid the double counting of the bottom friction in
+the 2-D barotropic momentum equations.
+As NEMO only updates the barotropic pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation,
+we need to remove the bottom friction induced by these two terms which has been included in the 3-D momentum trend
+and update it with the latest value.
+On the other hand, the bottom friction contributed by the other terms
+(\eg the advection term, viscosity term) has been included in the 3-D momentum equations and
+should not be added in the 2-D barotropic mode.
+The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the following:
+\[
+  % \label{eq:dynspg_ts_bfr1}
+  \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
+  \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
+\]
+\[
+  \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
+  \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
+\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
+\]
+where $\textbf{T}$ is the vertical integrated 3-D momentum trend.
+We assume the leap-frog time-stepping is used here.
+$\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
+$c_{b}$ is the friction coefficient.
+$\eta$ is the sea surface level calculated in the barotropic loops while $\eta^{'}$ is the sea surface level used in
+the 3-D baroclinic mode.
+$\textbf{u}_{b}$ is the bottom layer horizontal velocity.
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction with split-explicit time splitting
+% -------------------------------------------------------------------------------------------------------------
+\subsection[Bottom friction with split-explicit time splitting (\texttt{ln\_bfrimp})]
+{Bottom friction with split-explicit time splitting (\protect\np{ln\_bfrimp})}
+\label{subsec:ZDF_bfr_ts}
+When calculating the momentum trend due to bottom friction in \mdl{dynbfr},
+the bottom velocity at the before time step is used.
+This velocity includes both the baroclinic and barotropic components which is appropriate when
+using either the explicit or filtered surface pressure gradient algorithms
+(\key{dynspg\_exp} or \key{dynspg\_flt}).
+Extra attention is required, however, when using split-explicit time stepping (\key{dynspg\_ts}).
+In this case the free surface equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro},
+while the three dimensional prognostic variables are solved with the longer time step of \np{rn\_rdt} seconds.
+The trend in the barotropic momentum due to bottom friction appropriate to this method is that given by
+the selected parameterisation (\ie linear or non-linear bottom friction) computed with
+the evolving velocities at each barotropic timestep.
+In the case of non-linear bottom friction, we have elected to partially linearise the problem by
+keeping the coefficients fixed throughout the barotropic time-stepping to those computed in
+\mdl{zdfbfr} using the now timestep.
+This decision allows an efficient use of the $c_b^{\vect{U}}$ coefficients to:
+We recommend this option for shelf sea and coastal ocean applications. %, especially for split-explicit time splitting.
+This option can be invoked by setting \np{ln\_drgimp} to \forcode{.true.} in the \textit{namdrg} namelist.
+%This option requires \np{ln\_zdfexp} to be \forcode{.false.} in the \textit{namzdf} namelist.
+This implementation is performed in \mdl{dynzdf} where the following boundary conditions are set while solving the fully implicit diffusion step:
+At the top (below an ice shelf cavity):
+\[
+  % \label{eq:dynzdf_drg_top}
+  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{t}
+  = c_{t}^{\textbf{U}}\textbf{u}^{n+1}_{t}
+\]
+At the bottom (above the sea floor):
+\[
+  % \label{eq:dynzdf_drg_bot}
+  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{b}
+  = c_{b}^{\textbf{U}}\textbf{u}^{n+1}_{b}
+\]
+where $t$ and $b$ refers to top and bottom layers respectively.
+Superscript $n+1$ means the velocity used in the friction formula is to be calculated, so it is implicit.
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction with split-explicit free surface
+% -------------------------------------------------------------------------------------------------------------
+ \subsection[Bottom friction with split-explicit free surface]
+ {Bottom friction with split-explicit free surface}
+ \label{subsec:ZDF_drg_ts}
+With split-explicit free surface, the sub-stepping of barotropic equations needs the knowledge of top/bottom stresses. An obvious way to satisfy this is to take them as constant over the course of the barotropic integration and equal to the value used to update the baroclinic momentum trend. Provided \forcode{ln_drgimp = .false.} and a centred or \textit{leap-frog} like integration of barotropic equations is used (\ie \forcode{ln_bt_fw = .false.}, cf \autoref{subsec:DYN_spg_ts}), this does ensure that barotropic and baroclinic dynamics feel the same stresses during one leapfrog time step. However, if \forcode{ln_drgimp = .true.},  stresses depend on the \textit{after} value of the velocities which themselves depend on the barotropic iteration result. This cyclic dependency makes difficult obtaining consistent stresses in 2d and 3d dynamics. Part of this mismatch is then removed when setting the final barotropic component of 3d velocities to the time splitting estimate. This last step can be seen as a necessary evil but should be minimized since it interferes with the adjustment to the boundary conditions.
+The strategy to handle top/bottom stresses with split-explicit free surface in \NEMO is as follows:
 \begin{enumerate}
+\item On entry to \rou{dyn\_spg\_ts}, remove the contribution of the before barotropic velocity to
+  the bottom friction component of the vertically integrated momentum trend.
+  Note the same stability check that is carried out on the bottom friction coefficient in \mdl{dynbfr} has to
+  be applied here to ensure that the trend removed matches that which was added in \mdl{dynbfr}.
+\item At each barotropic step, compute the contribution of the current barotropic velocity to
+  the trend due to bottom friction.
+  Add this contribution to the vertically integrated momentum trend.
+  This contribution is handled implicitly which eliminates the need to impose a stability criteria on
+  the values of the bottom friction coefficient within the barotropic loop.
+\end{enumerate}
+Note that the use of an implicit formulation within the barotropic loop for the bottom friction trend means that
+any limiting of the bottom friction coefficient in \mdl{dynbfr} does not adversely affect the solution when
+using split-explicit time splitting.
+This is because the major contribution to bottom friction is likely to come from the barotropic component which
+uses the unrestricted value of the coefficient.
+However, if the limiting is thought to be having a major effect
+(a more likely prospect in coastal and shelf seas applications) then
+the fully implicit form of the bottom friction should be used (see \autoref{subsec:ZDF_bfr_imp})
+which can be selected by setting \np{ln\_bfrimp} $=$ \forcode{.true.}.
+Otherwise, the implicit formulation takes the form:
+\[
+  % \label{eq:zdfbfr_implicitts}
+  \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]
+\]
+where $\bar U$ is the barotropic velocity, $H_e$ is the full depth (including sea surface height),
+$c_b^u$ is the bottom friction coefficient as calculated in \rou{zdf\_bfr} and
+$RHS$ represents all the components to the vertically integrated momentum trend except for
+that due to bottom friction.
+% ================================================================
+% Tidal Mixing
+% ================================================================
+\section[Tidal mixing (\texttt{\textbf{key\_zdftmx}})]
+{Tidal mixing (\protect\key{zdftmx})}
+\label{sec:ZDF_tmx}
+%--------------------------------------------namzdf_tmx--------------------------------------------------
+\item To extend the stability of the barotropic sub-stepping, bottom stresses are refreshed at each sub-iteration. The baroclinic part of the flow entering the stresses is frozen at the initial time of the barotropic iteration. In case of non-linear friction, the drag coefficient is also constant.
+\item In case of an implicit drag, specific computations are performed in \mdl{dynzdf} which renders the overall scheme mixed explicit/implicit: the barotropic components of 3d velocities are removed before seeking for the implicit vertical diffusion result. Top/bottom stresses due to the barotropic components are explicitly accounted for thanks to the updated values of barotropic velocities. Then the implicit solution of 3d velocities is obtained. Lastly, the residual barotropic component is replaced by the time split estimate.
+\end{enumerate}
+Note that other strategies are possible, like considering vertical diffusion step in advance, \ie prior barotropic integration.
+% ================================================================
+% Internal wave-driven mixing
+% ================================================================
+\section[Internal wave-driven mixing (\forcode{ln_zdfiwm = .true.})]
+{Internal wave-driven mixing (\protect\np{ln\_zdfiwm}\forcode{ = .true.})}
+\label{subsec:ZDF_tmx_new}
+%--------------------------------------------namzdf_iwm------------------------------------------
+%
+%\nlst{namzdf_tmx}
+%--------------------------------------------------------------------------------------------------------------
+% -------------------------------------------------------------------------------------------------------------
+%        Bottom intensified tidal mixing
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom intensified tidal mixing}
+\label{subsec:ZDF_tmx_bottom}
+Options are defined through the  \ngn{namzdf\_tmx} namelist variables.
+The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by
+\citet{st-laurent.simmons.ea_GRL02} and first introduced in an OGCM by \citep{simmons.jayne.ea_OM04}.
+In this formulation an additional vertical diffusivity resulting from internal tide breaking,
+$A^{vT}_{tides}$ is expressed as a function of $E(x,y)$,
+the energy transfer from barotropic tides to baroclinic tides:
+\begin{equation}
+  \label{eq:Ktides}
+  A^{vT}_{tides} =  q \,\Gamma \,\frac{ E(x,y) \, F(z) }{ \rho \, N^2 }
+\end{equation}
+where $\Gamma$ is the mixing efficiency, $N$ the Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}),
+$\rho$ the density, $q$ the tidal dissipation efficiency, and $F(z)$ the vertical structure function.
+The mixing efficiency of turbulence is set by $\Gamma$ (\np{rn\_me} namelist parameter) and
+is usually taken to be the canonical value of $\Gamma = 0.2$ (Osborn 1980).
+The tidal dissipation efficiency is given by the parameter $q$ (\np{rn\_tfe} namelist parameter)
+represents the part of the internal wave energy flux $E(x, y)$ that is dissipated locally,
+with the remaining $1-q$ radiating away as low mode internal waves and
+contributing to the background internal wave field.
+A value of $q=1/3$ is typically used \citet{st-laurent.simmons.ea_GRL02}.
+The vertical structure function $F(z)$ models the distribution of the turbulent mixing in the vertical.
+It is implemented as a simple exponential decaying upward away from the bottom,
+with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter,
+with a typical value of $500\,m$) \citep{st-laurent.nash_DSR04},
+\[
+  % \label{eq:Fz}
+  F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
+\]
+and is normalized so that vertical integral over the water column is unity.
+The associated vertical viscosity is calculated from the vertical diffusivity assuming a Prandtl number of 1,
+\ie $A^{vm}_{tides}=A^{vT}_{tides}$.
+In the limit of $N \rightarrow 0$ (or becoming negative), the vertical diffusivity is capped at $300\,cm^2/s$ and
+impose a lower limit on $N^2$ of \np{rn\_n2min} usually set to $10^{-8} s^{-2}$.
+These bounds are usually rarely encountered.
+The internal wave energy map, $E(x, y)$ in \autoref{eq:Ktides}, is derived from a barotropic model of
+the tides utilizing a parameterization of the conversion of barotropic tidal energy into internal waves.
+The essential goal of the parameterization is to represent the momentum exchange between the barotropic tides and
+the unrepresented internal waves induced by the tidal flow over rough topography in a stratified ocean.
+In the current version of \NEMO, the map is built from the output of
+the barotropic global ocean tide model MOG2D-G \citep{carrere.lyard_GRL03}.
+This model provides the dissipation associated with internal wave energy for the M2 and K1 tides component
+(\autoref{fig:ZDF_M2_K1_tmx}).
+The S2 dissipation is simply approximated as being $1/4$ of the M2 one.
+The internal wave energy is thus : $E(x, y) = 1.25 E_{M2} + E_{K1}$.
+Its global mean value is $1.1$ TW,
+in agreement with independent estimates \citep{egbert.ray_N00, egbert.ray_JGR01}.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_ZDF_M2_K1_tmx}
+    \caption{
+      \protect\label{fig:ZDF_M2_K1_tmx}
+      (a) M2 and (b) K1 internal wave drag energy from \citet{carrere.lyard_GRL03} ($W/m^2$).
+    }
+  \end{center}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% -------------------------------------------------------------------------------------------------------------
+%        Indonesian area specific treatment
+% -------------------------------------------------------------------------------------------------------------
+\subsection[Indonesian area specific treatment (\texttt{ln\_zdftmx\_itf})]
+{Indonesian area specific treatment (\protect\np{ln\_zdftmx\_itf})}
+\label{subsec:ZDF_tmx_itf}
+When the Indonesian Through Flow (ITF) area is included in the model domain,
+a specific treatment of tidal induced mixing in this area can be used.
+It is activated through the namelist logical \np{ln\_tmx\_itf}, and the user must provide an input NetCDF file,
+\ifile{mask\_itf}, which contains a mask array defining the ITF area where the specific treatment is applied.
+When \np{ln\_tmx\_itf}\forcode{ = .true.}, the two key parameters $q$ and $F(z)$ are adjusted following
+the parameterisation developed by \citet{koch-larrouy.madec.ea_GRL07}:
+First, the Indonesian archipelago is a complex geographic region with a series of
+large, deep, semi-enclosed basins connected via numerous narrow straits.
+Once generated, internal tides remain confined within this semi-enclosed area and hardly radiate away.
+Therefore all the internal tides energy is consumed within this area.
+So it is assumed that $q = 1$, \ie all the energy generated is available for mixing.
+Note that for test purposed, the ITF tidal dissipation efficiency is a namelist parameter (\np{rn\_tfe\_itf}).
+A value of $1$ or close to is this recommended for this parameter.
+Second, the vertical structure function, $F(z)$, is no more associated with a bottom intensification of the mixing,
+but with a maximum of energy available within the thermocline.
+\citet{koch-larrouy.madec.ea_GRL07} have suggested that the vertical distribution of
+the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above.
+The resulting $F(z)$ is:
+\[
+  % \label{eq:Fz_itf}
+  F(i,j,k) \sim     \left\{
+    \begin{aligned}
+      \frac{q\,\Gamma E(i,j) } {\rho N \, \int N     dz}    \qquad \text{when $\partial_z N < 0$} \\
+      \frac{q\,\Gamma E(i,j) } {\rho     \, \int N^2 dz}    \qquad \text{when $\partial_z N > 0$}
+    \end{aligned}
+  \right.
+\]
+Averaged over the ITF area, the resulting tidal mixing coefficient is $1.5\,cm^2/s$,
+which agrees with the independent estimates inferred from observations.
+Introduced in a regional OGCM, the parameterization improves the water mass characteristics in
+the different Indonesian seas, suggesting that the horizontal and vertical distributions of
+the mixing are adequately prescribed \citep{koch-larrouy.madec.ea_GRL07, koch-larrouy.madec.ea_OD08*a, koch-larrouy.madec.ea_OD08*b}.
+Note also that such a parameterisation has a significant impact on the behaviour of
+global coupled GCMs \citep{koch-larrouy.lengaigne.ea_CD10}.
+% ================================================================
+% Internal wave-driven mixing
+% ================================================================
+\section[Internal wave-driven mixing (\texttt{\textbf{key\_zdftmx\_new}})]
+{Internal wave-driven mixing (\protect\key{zdftmx\_new})}
+\label{sec:ZDF_tmx_new}
+%--------------------------------------------namzdf_tmx_new------------------------------------------
+%
+%\nlst{namzdf_tmx_new}
+\nlst{namzdf_iwm}
 %--------------------------------------------------------------------------------------------------------------
 …
 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of
 the energy available for mixing.
 If the \np{ln\_mevar} namelist parameter is set to false, the mixing efficiency is taken as constant and
+If the \np{ln\_mevar} namelist parameter is set to \forcode{.false.}, the mixing efficiency is taken as constant and
 equal to 1/6 \citep{osborn_JPO80}.
 In the opposite (recommended) case, $R_f$ is instead a function of
 …
 In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary
 as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice.
+as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to \forcode{.true.}, a recommended choice.
 This parameterization of differential mixing, due to \cite{jackson.rehmann_JPO14},
 is implemented as in \cite{de-lavergne.madec.ea_JPO16}.
 …
 The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$,
 is constructed from three static maps of column-integrated internal wave energy dissipation,
 $E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures
+(de Lavergne et al., in prep):
+$E_{cri}(i,j)$, $E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures:
 \begin{align*}
   F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
   F_{pyc}(i,j,k) &\propto N^{n\_p}\\
+  F_{pyc}(i,j,k) &\propto N^{n_p}\\
   F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
 \end{align*}
 …
   h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
 \]
 The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)
+The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_iwm} namelist)
 controls the stratification-dependence of the pycnocline-intensified dissipation.
 It can take values of 1 (recommended) or 2.
+It can take values of $1$ (recommended) or $2$.
 Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of
 the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.
 …
 $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of
 the abyssal hill topography \citep{goff_JGR10} and the latitude.
+%
+% Jc: input files names ?
+% ================================================================
+% surface wave-induced mixing
+% ================================================================
+\section[Surface wave-induced mixing (\forcode{ln_zdfswm = .true.})]
+{Surface wave-induced mixing (\protect\np{ln\_zdfswm}\forcode{ = .true.})}
+\label{subsec:ZDF_swm}
+TBC ...
+% ================================================================
+% Adaptive-implicit vertical advection
+% ================================================================
+\section[Adaptive-implicit vertical advection (\forcode{ln_zad_Aimp = .true.})]
+{Adaptive-implicit vertical advection(\protect\np{ln\_zad\_Aimp}\forcode{ = .true.})}
+\label{subsec:ZDF_aimp}
+This refers to \citep{shchepetkin_OM15}.
+TBC ...
 % ================================================================

Note: See TracChangeset for help on using the changeset viewer.

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

Context Navigation

Changeset 11213

Legend:

NEMO/trunk/doc/latex/NEMO/main/bibliography.bib

NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

Download in other formats: