Changeset 10419 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_ZDF.tex

Timestamp:

2018-12-19T20:46:30+01:00 (5 years ago)

Author:

smasson

Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

Location:

NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex

Files:

• . (modified) (1 prop)
• NEMO (modified) (1 prop)
• NEMO/subfiles (modified) (1 prop)
• NEMO/subfiles/chap_ZDF.tex (modified) (53 diffs)

NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_ZDF.tex

@@ -1 @@
-\documentclass[../tex_main/NEMO_manual]{subfiles}
+\documentclass[../main/NEMO_manual]{subfiles}
+
 \begin{document}
 % ================================================================
@@ -6 +7 @@
 \chapter{Vertical Ocean Physics (ZDF)}
 \label{chap:ZDF}
+
 \minitoc
 
 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
 
-
 \newpage
-$\ $\newline    % force a new ligne
-
 
 % ================================================================
@@ -60 +59 @@
 It is recommended that this option is only used in process studies, not in basin scale simulations.
 Typical values used in this case are:
-\begin{align*} 
-A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}  \\
-A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
+\begin{align*}
+  A_u^{vm} = A_v^{vm} &= 1.2\ 10^{-4}~m^2.s^{-1}   \\
+  A^{vT} = A^{vS} &= 1.2\ 10^{-5}~m^2.s^{-1}
 \end{align*}
 
@@ -69 +68 @@
 that is $\sim10^{-6}~m^2.s^{-1}$ for momentum, $\sim10^{-7}~m^2.s^{-1}$ for temperature and
 $\sim10^{-9}~m^2.s^{-1}$ for salinity.
-
 
 % -------------------------------------------------------------------------------------------------------------
@@ -90 +88 @@
 ($i.e.$ the ratio of stratification to vertical shear).
 Following \citet{Pacanowski_Philander_JPO81}, the following formulation has been implemented:
-\begin{equation} \label{eq:zdfric}
-   \left\{      \begin{aligned}
-         A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
-         A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
-   \end{aligned}    \right.
-\end{equation}
+\[
+  % \label{eq:zdfric}
+  \left\{
+    \begin{aligned}
+      A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+a \; Ri \right)^n} + A_b^{vT}       \\
+      A^{vm} &= \frac{A^{vT}        }{\left( 1+ a \;Ri  \right)   } + A_b^{vm}
+    \end{aligned}
+  \right.
+\]
 where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number,
 $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 
@@ -111 +112 @@
 
 This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
-\begin{equation}
-h_{e} = Ek \frac {u^{*}} {f_{0}}
-\end{equation}
+\[
+  h_{e} = Ek \frac {u^{*}} {f_{0}}
+\]
 where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis parameter.
 
 In this similarity height relationship, the turbulent friction velocity:
-\begin{equation}
-u^{*} = \sqrt \frac {|\tau|} {\rho_o}
-\end{equation}
+\[
+  u^{*} = \sqrt \frac {|\tau|} {\rho_o}
+\]
 is computed from the wind stress vector $|\tau|$ and the reference density $ \rho_o$.
 The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
@@ -146 +147 @@
 The time evolution of $\bar{e}$ is the result of the production of $\bar{e}$ through vertical shear,
 its destruction through stratification, its vertical diffusion, and its dissipation of \citet{Kolmogorov1942} type:
-\begin{equation} \label{eq:zdftke_e}
-\frac{\partial \bar{e}}{\partial t} = 
-\frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
-                    +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
--K_\rho\,N^2
-+\frac{1}{e_3} \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
-            \;\frac{\partial \bar{e}}{\partial k}} \right]
-- c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
+\begin{equation}
+  \label{eq:zdftke_e}
+  \frac{\partial \bar{e}}{\partial t} =
+  \frac{K_m}{{e_3}^2 }\;\left[ {\left( {\frac{\partial u}{\partial k}} \right)^2
+      +\left( {\frac{\partial v}{\partial k}} \right)^2} \right]
+  -K_\rho\,N^2
+  +\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{A^{vm}}{e_3 }
+      \;\frac{\partial \bar{e}}{\partial k}} \right]
+  - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }
 \end{equation}
-\begin{equation} \label{eq:zdftke_kz}
-   \begin{split}
-         K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }     \\
-         K_\rho &= A^{vm} / P_{rt}
-   \end{split}
-\end{equation}
+\[
+  % \label{eq:zdftke_kz}
+  \begin{split}
+    K_m &= C_k\  l_k\  \sqrt {\bar{e}\; }    \\
+    K_\rho &= A^{vm} / P_{rt}
+  \end{split}
+\]
 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}), 
 $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, 
@@ -168 +171 @@
 They are set through namelist parameters \np{nn\_ediff} and \np{nn\_ediss}.
 $P_{rt}$ can be set to unity or, following \citet{Blanke1993}, be a function of the local Richardson number, $R_i$:
-\begin{align*} \label{eq:prt}
-P_{rt} = \begin{cases}
-                    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}  \\
-                    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}  \\
-                    \ \ 10 &      \text{if $\ 2 \leq R_i$} 
-            \end{cases}
+\begin{align*}
+  % \label{eq:prt}
+  P_{rt} =
+  \begin{cases}
+    \ \ \ 1 &      \text{if $\ R_i \leq 0.2$}   \\
+    5\,R_i &      \text{if $\ 0.2 \leq R_i \leq 2$}   \\
+    \ \ 10 &      \text{if $\ 2 \leq R_i$}
+  \end{cases}
 \end{align*}
 Options are defined through the  \ngn{namzdfy\_tke} namelist variables.
@@ -195 +200 @@
 
 \subsubsection{Turbulent length scale}
+
 For computational efficiency, the original formulation of the turbulent length scales proposed by
 \citet{Gaspar1990} has been simplified.
 Four formulations are proposed, the choice of which is controlled by the \np{nn\_mxl} namelist parameter.
 The first two are based on the following first order approximation \citep{Blanke1993}:
-\begin{equation} \label{eq:tke_mxl0_1}
-l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
+\begin{equation}
+  \label{eq:tke_mxl0_1}
+  l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N
 \end{equation}
 which is valid in a stable stratified region with constant values of the Brunt-Vais\"{a}l\"{a} frequency.
@@ -211 +218 @@
 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:
-\begin{equation} \label{eq:tke_mxl_constraint}
-\frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
-\qquad \text{with }\  l =  l_k = l_\epsilon
+\begin{equation}
+  \label{eq:tke_mxl_constraint}
+  \frac{1}{e_3 }\left| {\frac{\partial l}{\partial k}} \right| \leq 1
+  \qquad \text{with }\  l =  l_k = l_\epsilon
 \end{equation}
 \autoref{eq:tke_mxl_constraint} means that the vertical variations of the length scale cannot be larger than
@@ -227 +235 @@
 (and note that here we use numerical indexing):
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t] \begin{center}
-\includegraphics[width=1.00\textwidth]{Fig_mixing_length}
-\caption{ \protect\label{fig:mixing_length} 
-Illustration of the mixing length computation. }
-\end{center}  
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=1.00\textwidth]{Fig_mixing_length}
+    \caption{
+      \protect\label{fig:mixing_length}
+      Illustration of the mixing length computation.
+    }
+  \end{center}
 \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{equation} \label{eq:tke_mxl2}
-\begin{aligned}
- l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
+\[
+  % \label{eq:tke_mxl2}
+  \begin{aligned}
+    l_{up\ \ }^{(k)} &= \min \left(  l^{(k)} \ , \ l_{up}^{(k+1)} + e_{3t}^{(k)}\ \ \ \;  \right)
     \quad &\text{ from $k=1$ to $jpk$ }\ \\
- l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)   
+    l_{dwn}^{(k)} &= \min \left(  l^{(k)} \ , \ l_{dwn}^{(k-1)} + e_{3t}^{(k-1)}  \right)
     \quad &\text{ from $k=jpk$ to $1$ }\ \\
-\end{aligned}
-\end{equation}
+  \end{aligned}
+\]
 where $l^{(k)}$ is computed using \autoref{eq:tke_mxl0_1}, $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.
 
@@ -247 +259 @@
 $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the \np{nn\_mxl}\forcode{ = 3} case,
 the dissipation and mixing turbulent length scales are give as in \citet{Gaspar1990}:
-\begin{equation} \label{eq:tke_mxl_gaspar}
-\begin{aligned}
-& l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }    \\
-& l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right) 
-\end{aligned}
-\end{equation}
+\[
+  % \label{eq:tke_mxl_gaspar}
+  \begin{aligned}
+    & l_k          = \sqrt{\  l_{up} \ \ l_{dwn}\ }   \\
+    & l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)
+  \end{aligned}
+\]
 
 At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist parameter.
@@ -261 +274 @@
 $\bar{e}$ reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
 
-
 \subsubsection{Surface wave breaking parameterization}
 %-----------------------------------------------------------------------%
+
 Following \citet{Mellor_Blumberg_JPO04}, the TKE turbulence closure model has been modified to
 include the effect of surface wave breaking energetics.
@@ -272 +285 @@
 
 Following \citet{Craig_Banner_JPO94}, the boundary condition on surface TKE value is :
-\begin{equation}  \label{eq:ZDF_Esbc}
-\bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
+\begin{equation}
+  \label{eq:ZDF_Esbc}
+  \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
 \end{equation}
 where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality which depends on the ''wave age'',
 ranging from 57 for mature waves to 146 for younger waves \citep{Mellor_Blumberg_JPO04}. 
 The boundary condition on the turbulent length scale follows the Charnock's relation:
-\begin{equation} \label{eq:ZDF_Lsbc}
-l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
+\begin{equation}
+  \label{eq:ZDF_Lsbc}
+  l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
 \end{equation}
 where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
@@ -293 +308 @@
 surface $\bar{e}$ value.
 
-
 \subsubsection{Langmuir cells}
 %--------------------------------------%
+
 Langmuir circulations (LC) can be described as ordered large-scale vertical motions in
 the surface layer of the oceans.
@@ -313 +328 @@
 By making an analogy with the characteristic convective velocity scale ($e.g.$, \citet{D'Alessio_al_JPO98}),
 $P_{LC}$ is assumed to be : 
-\begin{equation}
+\[
 P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}
-\end{equation}
+\]
 where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
 With no information about the wave field, $w_{LC}$ is assumed to be proportional to 
@@ -322 +337 @@
   $u_s =  0.016 \,|U_{10m}|$.
   Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
-  $1.5~10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress}.
+  $1.5~10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress
+}.
 For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at
 a finite depth $H_{LC}$ (which is often close to the mixed layer depth),
 and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). 
 The resulting expression for $w_{LC}$ is :
-\begin{equation}
-w_{LC}  = \begin{cases}
-                   c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$}    \\
-                   0                             &      \text{otherwise} 
-                 \end{cases}
-\end{equation}
+\[
+  w_{LC}  =
+  \begin{cases}
+    c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$}    \\
+    0                             &      \text{otherwise}
+  \end{cases}
+\]
 where $c_{LC} = 0.15$ has been chosen by \citep{Axell_JGR02} as a good compromise to fit LES data.
 The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second.
@@ -341 +358 @@
 $H_{LC}$ is 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 
-\begin{equation}
+\[
 - \int_{-H_{LC}}^0 { N^2\;z  \;dz} = \frac{1}{2} u_s^2
-\end{equation}
-
+\]
 
 \subsubsection{Mixing just below the mixed layer}
@@ -362 +378 @@
 swell and waves is parameterized by \autoref{eq:ZDF_Esbc} the standard TKE surface boundary condition,
 plus a depth depend one given by:
-\begin{equation}  \label{eq:ZDF_Ehtau}
-S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau} 
+\begin{equation}
+  \label{eq:ZDF_Ehtau}
+  S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau}
 \end{equation}
 where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that
@@ -380 +397 @@
 They will be removed in the next release. 
 
-
-
 % from Burchard et al OM 2008 : 
 % the most critical process not reproduced by statistical turbulence models is the activity of 
@@ -390 +405 @@
 % (e.g. Mellor, 1989; Large et al., 1994; Meier, 2001; Axell, 2002; St. Laurent and Garrett, 2002).
 
-
-
 % -------------------------------------------------------------------------------------------------------------
 %        TKE discretization considerations
@@ -399 +412 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]   \begin{center}
-\includegraphics[width=1.00\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}  
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=1.00\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}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -419 +435 @@
 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}
+\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$,
@@ -443 +460 @@
 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}
+\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)$. 
@@ -468 +486 @@
 
 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}
+\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)
@@ -511 +530 @@
 $k$-$\omega$ \citep{Wilcox_1988} among others \citep{Umlauf_Burchard_JMS03,Kantha_Carniel_CSR05}). 
 The GLS scheme is given by the following set of equations:
-\begin{equation} \label{eq:zdfgls_e}
-\frac{\partial \bar{e}}{\partial t} = 
-\frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
-                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
--K_\rho \,N^2
-+\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
-- \epsilon
+\begin{equation}
+  \label{eq:zdfgls_e}
+  \frac{\partial \bar{e}}{\partial t} =
+  \frac{K_m}{\sigma_e e_3 }\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
+      +\left( \frac{\partial v}{\partial k} \right)^2} \right]
+  -K_\rho \,N^2
+  +\frac{1}{e_3}\,\frac{\partial}{\partial k} \left[ \frac{K_m}{e_3}\,\frac{\partial \bar{e}}{\partial k} \right]
+  - \epsilon
 \end{equation}
 
-\begin{equation} \label{eq:zdfgls_psi}
-   \begin{split}
-\frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
-\frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
-                                                                   +\left( \frac{\partial v}{\partial k} \right)^2} \right]
-- C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
-&+\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
-                                  \;\frac{\partial \psi}{\partial k}} \right]\;
-   \end{split}
-\end{equation}
-
-\begin{equation} \label{eq:zdfgls_kz}
-   \begin{split}
-         K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
-         K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
-   \end{split}
-\end{equation}
-
-\begin{equation} \label{eq:zdfgls_eps}
-{\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
-\end{equation}
+\[
+  % \label{eq:zdfgls_psi}
+  \begin{split}
+    \frac{\partial \psi}{\partial t} =& \frac{\psi}{\bar{e}} \left\{
+      \frac{C_1\,K_m}{\sigma_{\psi} {e_3}}\;\left[ {\left( \frac{\partial u}{\partial k} \right)^2
+          +\left( \frac{\partial v}{\partial k} \right)^2} \right]
+      - C_3 \,K_\rho\,N^2   - C_2 \,\epsilon \,Fw   \right\}             \\
+    &+\frac{1}{e_3}  \;\frac{\partial }{\partial k}\left[ {\frac{K_m}{e_3 }
+        \;\frac{\partial \psi}{\partial k}} \right]\;
+  \end{split}
+\]
+
+\[
+  % \label{eq:zdfgls_kz}
+  \begin{split}
+    K_m    &= C_{\mu} \ \sqrt {\bar{e}} \ l         \\
+    K_\rho &= C_{\mu'}\ \sqrt {\bar{e}} \ l
+  \end{split}
+\]
+
+\[
+  % \label{eq:zdfgls_eps}
+  {\epsilon} = C_{0\mu} \,\frac{\bar {e}^{3/2}}{l} \;
+\]
 where $N$ is the local Brunt-Vais\"{a}l\"{a} frequency (see \autoref{subsec:TRA_bn2}) and
 $\epsilon$ the dissipation rate. 
@@ -549 +572 @@
 
 %--------------------------------------------------TABLE--------------------------------------------------
-\begin{table}[htbp]  \begin{center}
-%\begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
-\begin{tabular}{ccccc}
-                         &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\  
-%                        & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\  
-\hline  \hline 
-\np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\  
-\hline 
-$( p , n , m )$          &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
-$\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
-$\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
-$C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
-$C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
-$C_3$              &      1.           &     1.              &      1.                &       1.           \\
-$F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
-\hline
-\hline
-\end{tabular}
-\caption{   \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}.}
-\end{center}   \end{table}
+\begin{table}[htbp]
+  \begin{center}
+    % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}cp{70pt}c}
+    \begin{tabular}{ccccc}
+      &   $k-kl$   & $k-\epsilon$ & $k-\omega$ &   generic   \\
+      % & \citep{Mellor_Yamada_1982} &  \citep{Rodi_1987}       & \citep{Wilcox_1988} &                 \\
+      \hline
+      \hline
+      \np{nn\_clo}     & \textbf{0} &   \textbf{1}  &   \textbf{2}   &    \textbf{3}   \\
+      \hline
+      $( p , n , m )$          &   ( 0 , 1 , 1 )   & ( 3 , 1.5 , -1 )   & ( -1 , 0.5 , -1 )    &  ( 2 , 1 , -0.67 )  \\
+      $\sigma_k$      &    2.44         &     1.              &      2.                &      0.8          \\
+      $\sigma_\psi$  &    2.44         &     1.3            &      2.                 &       1.07       \\
+      $C_1$              &      0.9         &     1.44          &      0.555          &       1.           \\
+      $C_2$              &      0.5         &     1.92          &      0.833          &       1.22       \\
+      $C_3$              &      1.           &     1.              &      1.                &       1.           \\
+      $F_{wall}$        &      Yes        &       --             &     --                  &      --          \\
+      \hline
+      \hline
+    \end{tabular}
+    \caption{
+      \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}.
+    }
+  \end{center}
+\end{table}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -646 +674 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!htb]    \begin{center}
-\includegraphics[width=0.90\textwidth]{Fig_npc}
-\caption{  \protect\label{fig:npc} 
-  Example of an unstable density profile treated by the non penetrative convective adjustment algorithm.
-  $1^{st}$ step: the initial profile is checked from the surface to the bottom.
-  It is found to be unstable between levels 3 and 4.
-  They are mixed.
-  The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed.
-  The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed.
-  The $1^{st}$ step ends since the density profile is then stable below the level 3.
-  $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step:
-  levels 2 to 5 are mixed.
-  The new density profile is checked.
-  It is found stable: end of algorithm.}
-\end{center}   \end{figure}
+\begin{figure}[!htb]
+  \begin{center}
+    \includegraphics[width=0.90\textwidth]{Fig_npc}
+    \caption{
+      \protect\label{fig:npc}
+      Example of an unstable density profile treated by the non penetrative convective adjustment algorithm.
+      $1^{st}$ step: the initial profile is checked from the surface to the bottom.
+      It is found to be unstable between levels 3 and 4.
+      They are mixed.
+      The resulting $\rho$ is still larger than $\rho$(5): levels 3 to 5 are mixed.
+      The resulting $\rho$ is still larger than $\rho$(6): levels 3 to 6 are mixed.
+      The $1^{st}$ step ends since the density profile is then stable below the level 3.
+      $2^{nd}$ step: the new $\rho$ profile is checked following the same procedure as in $1^{st}$ step:
+      levels 2 to 5 are mixed.
+      The new density profile is checked.
+      It is found stable: end of algorithm.
+    }
+  \end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -781 +813 @@
 
 Diapycnal mixing of S and T are described by diapycnal diffusion coefficients 
-\begin{align*} % \label{eq:zdfddm_Kz}
-    &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT}  \\
-    &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
+\begin{align*}
+  % \label{eq:zdfddm_Kz}
+  &A^{vT} = A_o^{vT}+A_f^{vT}+A_d^{vT} \\
+  &A^{vS} = A_o^{vS}+A_f^{vS}+A_d^{vS}
 \end{align*}
 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection,
@@ -792 +825 @@
 To represent mixing of $S$ and $T$ by salt fingering, we adopt the diapycnal diffusivities suggested by Schmitt
 (1981):
-\begin{align} \label{eq:zdfddm_f}
-A_f^{vS} &=    \begin{cases}
-   \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
-   0                              &\text{otherwise} 
-            \end{cases}   
-\\           \label{eq:zdfddm_f_T}
-A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
+\begin{align}
+  \label{eq:zdfddm_f}
+  A_f^{vS} &=
+             \begin{cases}
+               \frac{A^{\ast v}}{1+(R_\rho / R_c)^n   } &\text{if  $R_\rho > 1$ and $N^2>0$ } \\
+               0                              &\text{otherwise}
+             \end{cases}
+  \\         \label{eq:zdfddm_f_T}
+  A_f^{vT} &= 0.7 \ A_f^{vS} / R_\rho 
 \end{align}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]   \begin{center}
-\includegraphics[width=0.99\textwidth]{Fig_zdfddm}
-\caption{  \protect\label{fig:zdfddm}
-  From \citet{Merryfield1999} :
-  (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering.
-  Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
-  (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in regions of diffusive convection.
-  Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation.
-  The latter is not implemented in \NEMO. }
-\end{center}    \end{figure}
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=0.99\textwidth]{Fig_zdfddm}
+    \caption{
+      \protect\label{fig:zdfddm}
+      From \citet{Merryfield1999} :
+      (a) Diapycnal diffusivities $A_f^{vT}$ and $A_f^{vS}$ for temperature and salt in regions of salt fingering.
+      Heavy curves denote $A^{\ast v} = 10^{-3}~m^2.s^{-1}$ and thin curves $A^{\ast v} = 10^{-4}~m^2.s^{-1}$;
+      (b) diapycnal diffusivities $A_d^{vT}$ and $A_d^{vS}$ for temperature and salt in regions of
+      diffusive convection.
+      Heavy curves denote the Federov parameterisation and thin curves the Kelley parameterisation.
+      The latter is not implemented in \NEMO.
+    }
+  \end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
@@ -820 +860 @@
 To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by
 Federov (1988) is used: 
-\begin{align}  \label{eq:zdfddm_d}
-A_d^{vT} &=    \begin{cases}
-   1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
-                           &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
-   0                       &\text{otherwise} 
-            \end{cases}   
-\\          \label{eq:zdfddm_d_S}
-A_d^{vS} &=    \begin{cases}
-   A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right)
-                           &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
-   A_d^{vT} \ 0.15 \ R_\rho
-                           &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
-   0                       &\text{otherwise} 
-            \end{cases}   
+\begin{align}
+  % \label{eq:zdfddm_d}
+  A_d^{vT} &=
+             \begin{cases}
+               1.3635 \, \exp{\left( 4.6\, \exp{ \left[  -0.54\,( R_{\rho}^{-1} - 1 )  \right] }    \right)}
+               &\text{if  $0<R_\rho < 1$ and $N^2>0$ } \\
+               0                       &\text{otherwise}
+             \end{cases}
+                                       \nonumber \\
+  \label{eq:zdfddm_d_S}
+  A_d^{vS} &=
+             \begin{cases}
+               A_d^{vT}\ \left( 1.85\,R_{\rho} - 0.85 \right) &\text{if  $0.5 \leq R_\rho<1$ and $N^2>0$ } \\
+               A_d^{vT} \ 0.15 \ R_\rho               &\text{if  $\ \ 0 < R_\rho<0.5$ and $N^2>0$ } \\
+               0                       &\text{otherwise}
+             \end{cases}
 \end{align}
 
@@ -863 +905 @@
 a condition on the vertical diffusive flux.
 For the bottom boundary layer, one has:
-\begin{equation} \label{eq:zdfbfr_flux}
-A^{vm} \left( \partial {\textbf U}_h / \partial z \right) = {{\cal F}}_h^{\textbf U}
-\end{equation}
+\[
+  % \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).
@@ -878 +921 @@
 bottom model layer.
 To illustrate this, consider the equation for $u$ at $k$, the last ocean level:
-\begin{equation} \label{eq:zdfbfr_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}}
+\begin{equation}
+  \label{eq:zdfbfr_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}
 If the bottom layer thickness is 200~m, the Ekman transport will be distributed over that depth.
@@ -897 +941 @@
 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:
-\begin{equation} \label{eq:zdfbfr_bdef}
-\frac{\partial {\textbf U_h}}{\partial t} = 
+\begin{equation}
+  \label{eq:zdfbfr_bdef}
+  \frac{\partial {\textbf U_h}}{\partial t} =
   - \frac{{\cal F}^{\textbf U}_{h}}{e_{3u}} = \frac{c_b^{\textbf U}}{e_{3u}} \;{\textbf U}_h^b
 \end{equation}
@@ -911 +956 @@
 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 (i.e. the velocity of the last model level):
-\begin{equation} \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
-\end{equation}
+\[
+  % \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}$.
 This coefficient is generally estimated by setting a typical decay time $\tau$ in the deep ocean, 
@@ -927 +973 @@
 
 For the linear friction case the coefficients defined in the general expression \autoref{eq:zdfbfr_bdef} are: 
-\begin{equation} \label{eq:zdfbfr_linbfr_b}
-\begin{split}
- c_b^u &= - r\\
- c_b^v &= - r\\
-\end{split}
-\end{equation}
+\[
+  % \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
@@ -950 +997 @@
 
 The non-linear bottom friction parameterisation assumes that the bottom friction is quadratic: 
-\begin{equation} \label{eq:zdfbfr_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 
-\end{equation}
+\[
+  % \label{eq:zdfbfr_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,
 internal waves breaking and other short time scale currents.
@@ -965 +1013 @@
 the bottom friction to the general momentum trend in \mdl{dynbfr}.
 For the non-linear friction case the terms computed in \mdl{zdfbfr} are:
-\begin{equation} \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}
-\end{equation}
+\[
+  % \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:
@@ -991 +1040 @@
 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:
-\begin{equation}
-C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
-\end{equation}
+\[
+  C_D = \left ( {\kappa \over {\rm 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
@@ -1001 +1050 @@
 the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional namelist parameter:
 \np{rn\_bfri2\_max}, i.e.:
-\begin{equation}
-rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
-\end{equation}
+\[
+  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
@@ -1018 +1067 @@
 bottom friction does not induce numerical instability.
 For the purposes of stability analysis, an approximation to \autoref{eq:zdfbfr_flux2} is:
-\begin{equation} \label{eq:Eqn_bfrstab}
-\begin{split}
- \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
-               &= -\frac{ru}{e_{3u}}\;2\rdt\\
-\end{split}
+\begin{equation}
+  \label{eq:Eqn_bfrstab}
+  \begin{split}
+    \Delta u &= -\frac{{{\cal F}_h}^u}{e_{3u}}\;2 \rdt    \\
+    &= -\frac{ru}{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:
-\begin{equation}
- |\Delta u| < \;|u| 
-\end{equation}
+\[
+  |\Delta u| < \;|u| 
+\]
 \noindent which, using \autoref{eq:Eqn_bfrstab}, gives:
-\begin{equation}
-r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
-\end{equation}
+\[
+  r\frac{2\rdt}{e_{3u}} < 1 \qquad  \Rightarrow \qquad r < \frac{e_{3u}}{2\rdt}\\
+\]
 This same inequality can also be derived in the non-linear bottom friction case if
 a velocity of 1 m.s$^{-1}$ is assumed.
 Alternatively, this criterion can be rearranged to suggest a minimum bottom box thickness to ensure stability:
-\begin{equation}
-e_{3u} > 2\;r\;\rdt
-\end{equation}
+\[
+  e_{3u} > 2\;r\;\rdt
+\]
 \noindent which it may be necessary to impose if partial steps are being used.
 For example, if $|u| = 1$ m.s$^{-1}$, $rdt = 1800$ s, $r = 10^{-3}$ then $e_{3u}$ should be greater than 3.6 m.
@@ -1070 +1120 @@
 the bottom boundary condition is implemented implicitly.
 
-\begin{equation} \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}}
-\end{equation}
+\[
+  % \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.
@@ -1089 +1140 @@
 The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the following:
 
-\begin{equation} \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)
-\end{equation}
-\begin{equation} \label{eq:dynspg_ts_bfr2}
-\frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
-\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
-2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
-\end{equation}
+\[
+  % \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)-
+  2\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.
@@ -1106 +1158 @@
 the 3-D baroclinic mode.
 $\textbf{u}_{b}$ is the bottom layer horizontal velocity.
-
-
-
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1157 +1206 @@
 
 Otherwise, the implicit formulation takes the form:
-\begin{equation} \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 ]  
-\end{equation}
+\[
+  % \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.
-
-
-
 
 % ================================================================
@@ -1192 +1239 @@
 $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 }
+\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}),
@@ -1209 +1257 @@
 with a vertical scale of $h_o$ (\np{rn\_htmx} namelist parameter,
 with a typical value of $500\,m$) \citep{St_Laurent_Nash_DSR04}, 
-\begin{equation} \label{eq:Fz}
-F(i,j,k) = \frac{ e^{ -\frac{H+z}{h_o} } }{ h_o \left( 1- e^{ -\frac{H}{h_o} } \right) }
-\end{equation}
+\[
+  % \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. 
 
@@ -1234 +1283 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]   \begin{center}
-\includegraphics[width=0.90\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}
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=0.90\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}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  
@@ -1268 +1321 @@
 the energy dissipation proportional to $N^2$ below the core of the thermocline and to $N$ above. 
 The resulting $F(z)$ is:
-\begin{equation} \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.
-\end{equation}
+\[
+  % \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$, 
@@ -1283 +1339 @@
 global coupled GCMs \citep{Koch-Larrouy_al_CD10}.
 
-
 % ================================================================
 % Internal wave-driven mixing
@@ -1299 +1354 @@
 A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,
 and the resulting diffusivity is obtained as 
-\begin{equation} \label{eq:Kwave}
-A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
-\end{equation}
+\[
+  % \label{eq:Kwave}
+  A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
+\]
 where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution of
 the energy available for mixing.
@@ -1323 +1379 @@
 (de Lavergne et al., in prep):
 \begin{align*}
-F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
-F_{pyc}(i,j,k) &\propto N^{n\_p}\\
-F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
+  F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
+  F_{pyc}(i,j,k) &\propto N^{n\_p}\\
+  F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
 \end{align*} 
 In the above formula, $h_{ab}$ denotes the height above bottom,
 $h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
-\begin{equation*}
-h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
-\end{equation*}
+\[
+  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)
 controls the stratification-dependence of the pycnocline-intensified dissipation.
@@ -1343 +1399 @@
 % ================================================================
 
-
+\biblio
 
 \end{document}