Changeset 11422

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

r10075	r11422
22	22	cn_exp = "Pacific" ! experience name
23	23	nn_it000 = 1 ! first time step
24		nn_itend = 75 ! last time step (std 5475)
	24	nn_itend = 80 ! last time step (std 5475)
25	25	nn_istate = 0 ! output the initial state (1) or not (0)
26	26	/
…	…
30	30	ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time
31	31	!
32		rn_rdt = 5760. ! time step for the dynamics and tracer
	32	rn_rdt = 5400. ! time step for the dynamics and tracer
33	33	/
34	34	!-----------------------------------------------------------------------

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

-                      r10072
+                      r11422
    cn_exp      = "Nordic"  !  experience name
    nn_it000    =    1      !  first time step
    nn_itend    =  300      !  last  time step
+   nn_itend    =  320      !  last  time step
    nn_istate   =      0    !  output the initial state (1) or not (0)
    ln_clobber  = .true.    !  clobber (overwrite) an existing file
 …
    ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time
+   !
    rn_rdt      = 1440.     !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_rdt      = 1350.     !  time step for the dynamics (and tracer if nn_acc=0)
+/
 !-----------------------------------------------------------------------

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

-                      r10072
+                      r11422
    cn_exp      = "Nordic"  !  experience name
    nn_it000    =    1      !  first time step
    nn_itend    =  900      !  last  time step
+   nn_itend    =  960      !  last  time step
    nn_istate   =      0    !  output the initial state (1) or not (0)
    ln_clobber  = .true.    !  clobber (overwrite) an existing file
 …
    ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time
+   !
    rn_rdt      = 480.     !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_rdt      = 450.     !  time step for the dynamics (and tracer if nn_acc=0)
+/
 !-----------------------------------------------------------------------

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

r10075	r11422
22	22	cn_exp = "ORCA2" ! experience name
23	23	nn_it000 = 1 ! first time step
24		nn_itend = 75 ! last time step (std 5475)
	24	nn_itend = 80 ! last time step (std 5475)
25	25	nn_istate = 0 ! output the initial state (1) or not (0)
26	26	/
…	…
30	30	ln_linssh = .false. ! =T linear free surface ==>> model level are fixed in time
31	31	!
32		rn_rdt = 5760. ! time step for the dynamics and tracer
	32	rn_rdt = 5400. ! time step for the dynamics and tracer
33	33	/
34	34	!-----------------------------------------------------------------------

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

-                      r10190
+                      r11422
    cn_exp      =  "ORCA2"  !  experience name
    nn_it000    =       1   !  first time step
    nn_itend    =    5475   !  last  time step (std 5475)
+   nn_itend    =    5840   !  last  time step (std 5475)
    nn_istate   =       0   !  output the initial state (1) or not (0)
+/
 …
 &namdom        !   time and space domain
 !-----------------------------------------------------------------------
    rn_rdt      = 5760.     !  time step for the dynamics and tracer
+   rn_rdt      = 5400.     !  time step for the dynamics and tracer
+/
 !-----------------------------------------------------------------------
 …
 &namsbc        !   Surface Boundary Condition manager                   (default: NO selection)
 !-----------------------------------------------------------------------
    nn_fsbc     = 3         !  frequency of SBC module call
+   nn_fsbc     = 4         !  frequency of SBC module call
                            !     (also = the frequency of sea-ice & iceberg model call)
                      ! Type of air-sea fluxes

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

r10911	r11422
102	102	&namdia ! Diagnostics
103	103	!------------------------------------------------------------------------------
	104	ln_icediachk = .true. ! check online the heat, mass & salt budgets (T) or not (F)
104	105	/

NEMO/branches/2019/fix_vvl_ticket1791/doc

Property svn:ignore deleted

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex

Property svn:ignore

-                      old
+                      new
+*.aux
+*.bbl
+*.blg
+*.dvi
+*.fdb*
+*.fls
+*.idx
+*.ilg
+*.ind
+*.log
+*.maf
+*.mtc*
+*.out
+*.pdf
+*.toc
+_minted-*
+figures

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/HTML_latex2html.sh

-                      r10146
+                      r11422
 #!/bin/bash
 ./inc/clean.sh
 ./inc/build.sh
+#./inc/clean.sh
+#./inc/build.sh
+cd tex_main
+sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \
+   ../tex_sub/*.tex
+sed -i    's#\\subfile{#\\include{#g' \
+   NEMO_manual.tex
+latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML   \
+            $*                                                                         \
+   NEMO_manual
+sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \
+   ../tex_sub/*.tex
+sed -i    's#\\include{#\\subfile{#g' \
+   NEMO_manual.tex
+sed -i -e 's#utf8#latin1#'                                 \
+       -e 's#\[outputdir=../build\]{minted}#\[\]{minted}#' \
+       -e '/graphicspath/ s#{../#{../../#g'                \
+          global/packages.tex
+cd ./NEMO/main
+sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s#^#%#'   ../subfiles/*.tex
+sed -i    's#\\subfile{#\\input{#'                                         chapters.tex appendices.tex
+#latex2html                                   -noimages -local_icons -no_footnode -split 4 -link 2 -dir ../html_LaTeX2HTML $* NEMO_manual
+latex2html -debug -noreuse -init_file ../../l2hconf.pm -local_icons                               -dir ../build/html               NEMO_manual
+sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s#^%##'   ../subfiles/*.tex
+sed -i    's#\\input{#\\subfile{#'                                         chapters.tex appendices.tex
 cd -
+sed -i -e 's#latin1#utf8#'                                 \
+       -e 's#\[\]{minted}#\[outputdir=../build\]{minted}#' \
+       -e '/graphicspath/ s#{../../#{../#g'                \
+     global/packages.tex
 exit 0

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO

Property svn:ignore deleted

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/main

Property svn:ignore

-                      old
+                      new
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
 *.ilg
 *.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*

Property svn:externals set to
^/utils/dev/bibtool.rsc bibtool.rsc

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_A.tex

-                      r10442
+                      r11422
 \minitoc
+\vfill
+\begin{figure}[b]
+\subsubsection*{Changes record}
+\begin{tabular}{l||l|m{0.65\linewidth}}
+    Release   & Author        & Modifications \\
+    {\em 4.0} & {\em Mike Bell} & {\em review}  \\
+    {\em 3.x} & {\em Gurvan Madec} & {\em original}  \\
+\end{tabular}
+\end{figure}
 \newpage
 …
 \begin{equation}
   \label{apdx:A_s_slope}
   \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s
+  \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s
   \quad \text{and} \quad
+  \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s
+\end{equation}
+The chain rule to establish the model equations in the curvilinear $s-$coordinate system is:
+  \sigma_2 =\frac{1}{e_2 } \; \left. {\frac{\partial z}{\partial j}} \right|_s .
+\end{equation}
+The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as
+functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of
+these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms:
+\begin{equation}
+  \label{apdx:A_s_infin_changes}
+  \begin{aligned}
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t}
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t}
+                + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t}
+                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t}
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t}
+                + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t}
+                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} .
+  \end{aligned}
+\end{equation}
+Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that
+\begin{equation}
+  \label{apdx:A_s_chain_rule}
+      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t}  =
+      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t}
+    + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \;
+      \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .
+\end{equation}
+The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces,
+(\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to
+$s$ and $j, t$ held constant
+\begin{equation}
+\label{apdx:a_delta_s}
+\delta s|_{j,t} =
+         \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t}
+       + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} .
+\end{equation}
+Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using
+(\autoref{apdx:A_s_slope}) we obtain
+\begin{equation}
+\left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =
+         -  \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \;
+            \left. \frac{ \partial s }{\partial z} \right|_{i,j,t}
+    = - \frac{e_1 }{e_3 }\sigma_1  .
+\label{apdx:a_ds_di_z}
+\end{equation}
+Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived
+by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider
+changes in which $i , j$ and $s$ are constant. This shows that
+\begin{equation}
+\label{apdx:A_w_in_s}
+w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =
+- \left. \frac{ \partial z }{\partial s} \right|_{i,j,t}
+  \left. \frac{ \partial s }{\partial t} \right|_{i,j,z}
+  = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} .
+\end{equation}
+In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is
+usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish
+the model equations in the curvilinear $s-$coordinate system are:
 \begin{equation}
   \label{apdx:A_s_chain_rule}
   \begin{aligned}
     &\left. {\frac{\partial \bullet }{\partial t}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial t}} \right|_s
     -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial t} \\
+    \left. {\frac{\partial \bullet }{\partial t}} \right|_s
+    + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\
     &\left. {\frac{\partial \bullet }{\partial i}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial i}} \right|_s
     -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial i}=
     \left. {\frac{\partial \bullet }{\partial i}} \right|_s
     -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} \\
+    +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}=
+    \left. {\frac{\partial \bullet }{\partial i}} \right|_s
+    -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\
     &\left. {\frac{\partial \bullet }{\partial j}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial j}} \right|_s
     - \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=
     \left. {\frac{\partial \bullet }{\partial j}} \right|_s
     - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} \\
     &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s}
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
+    + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
+    - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\
+    &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} .
   \end{aligned}
 \end{equation}
-In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$,
-the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate:
-\begin{equation}
-  \label{apdx:A_w_in_s}
-  w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s
-  = \frac{\partial z}{\partial s} \; \frac{\partial s}{\partial t}
-  = e_3 \, \frac{\partial s}{\partial t}
-\end{equation}
 % ================================================================
 …
+    {
     \begin{array}{*{20}l}
       \nabla \cdot {\rm {\bf U}}
+      \nabla \cdot {\mathrm {\mathbf U}}
       &= \frac{1}{e_1 \,e_2 }  \left[ \left. {\frac{\partial (e_2 \,u)}{\partial i}} \right|_z
         +\left. {\frac{\partial(e_1 \,v)}{\partial j}} \right|_z  \right]
 …
       $, it becomes:}
+    %
       \nabla \cdot {\rm {\bf U}}
+      \nabla \cdot {\mathrm {\mathbf U}}
       & = \frac{1}{e_1 \,e_2 \,e_3 }  \left[
         \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
 …
 \end{subequations}
+Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
+Introducing the dia-surface velocity component,
+$\omega $, defined as the volume flux across the moving $s$-surfaces per unit horizontal area:
+Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
+Using the first form of (\autoref{apdx:A_s_infin_changes})
+and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$,
+one can show that the vertical velocity, $w_p$ of a point
+moving with the horizontal velocity of the fluid along an $s$ surface is given by
+\begin{equation}
+\label{apdx:A_w_p}
+\begin{split}
+w_p  = & \left. \frac{ \partial z }{\partial t} \right|_s
+     + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s
+     + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\
+     = & w_s + u \sigma_1 + v \sigma_2 .
+\end{split}
+\end{equation}
+ The vertical velocity across this surface is denoted by
 \begin{equation}
   \label{apdx:A_w_s}
+  \omega  = w - w_s - \sigma_1 \,u - \sigma_2 \,v    \\
+\end{equation}
+with $w_s$ given by \autoref{apdx:A_w_in_s},
+we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
+\begin{subequations}
+  \begin{align*}
+    {
+    \begin{array}{*{20}l}
+      \nabla \cdot {\rm {\bf U}}
+      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
+  \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  .
+\end{equation}
+Hence
+\begin{equation}
+\frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] =
+\frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] =
+   \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s}
+ + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] =
+   \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s
+\end{equation}
+Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain
+our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
+\begin{equation}
+      \nabla \cdot {\mathrm {\mathbf U}} =
+         \frac{1}{e_1 \,e_2 \,e_3 }    \left[
         \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
         +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
         + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
+        + \frac{1}{e_3 } \frac{\partial w_s       }{\partial s} \\ \\
+      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
+        \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
+        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
+        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
+        + \frac{1}{e_3 } \frac{\partial}{\partial s}  \left(  e_3 \; \frac{\partial s}{\partial t}   \right) \\ \\
+      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
+        \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
+        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
+        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
+        + \frac{\partial}{\partial s} \frac{\partial s}{\partial t}
+        + \frac{1}{e_3 } \frac{\partial s}{\partial t} \frac{\partial e_3}{\partial s} \\ \\
+      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
+        \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
+        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
+        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
+        + \frac{1}{e_3 } \frac{\partial e_3}{\partial t}
+    \end{array}
+        }
+  \end{align*}
+\end{subequations}
+        + \frac{1}{e_3 } \left. \frac{\partial e_3}{\partial t} \right|_s .
+\end{equation}
 As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is:
 …
     {\left. {\frac{\partial (e_2 \,e_3 \,u)}{\partial i}} \right|_s
       +  \left. {\frac{\partial (e_1 \,e_3 \,v)}{\partial j}} \right|_s } \right]
   +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0
 \end{equation}
 A additional term has appeared that take into account
+  +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 .
+\end{equation}
+An additional term has appeared that takes into account
 the contribution of the time variation of the vertical coordinate to the volume budget.
 …
         + w \;\frac{\partial u}{\partial z} \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
+        - \left. \zeta \right|_z v
+        +  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z
+        -  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z
         -\left.{ \frac{\partial (e_1 \,u)}{\partial j} }\right|_z } \right] \; v
         +  \frac{1}{2e_1} \left.{ \frac{\partial (u^2+v^2)}{\partial i} } \right|_z
 …
         } \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
         + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\
       &\qquad \qquad \qquad \quad
         + \frac{w}{e_3 } \;\frac{\partial u}{\partial s}
         - \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s}
+        + \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s}
         - \frac{\sigma_2 }{e_3 }\frac{\partial u}{\partial s}}   \right]\;v
         - \frac{\sigma_1 }{2e_3 }\frac{\partial (u^2+v^2)}{\partial s} \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
         + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\
       &\qquad \qquad \qquad \quad
 …
         - \sigma_1 u\frac{\partial u}{\partial s} - \sigma_1 v\frac{\partial v}{\partial s}} \right] \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
         + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
         + \frac{1}{e_3} \left[  w - \sigma_2 v - \sigma_1 u  \right]
         \; \frac{\partial u}{\partial s}   \\
+        \; \frac{\partial u}{\partial s} .  \\
+        %
       \intertext{Introducing $\omega$, the dia-a-surface velocity given by (\autoref{apdx:A_w_s}) }
+      \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) }
+      %
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
         + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
         + \frac{1}{e_3 } \left( \omega - w_s \right) \frac{\partial u}{\partial s}   \\
+        + \frac{1}{e_3 } \left( \omega + w_s \right) \frac{\partial u}{\partial s}   \\
     \end{array}
+    }
 …
+  {
     \begin{array}{*{20}l}
       w_s  \;\frac{\partial u}{\partial s}
       = \frac{\partial s}{\partial t} \;  \frac{\partial u }{\partial s}
       = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \quad ,
+      \frac{w_s}{e_3}  \;\frac{\partial u}{\partial s}
+      = - \left. \frac{\partial s}{\partial t} \right|_z \;  \frac{\partial u }{\partial s}
+      = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \ .
     \end{array}
+  }
 \]
 leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,
+This leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,
 \ie the total $s-$coordinate time derivative :
 \begin{align}
 …
   \left. \frac{D u}{D t} \right|_s
   = \left. {\frac{\partial u }{\partial t}} \right|_s
   + \left. \zeta \right|_s \;v
+  - \left. \zeta \right|_s \;v
   + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
   + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s}
+  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} .
 \end{align}
 Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in
 …
                                          + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
                                       &&- \frac{v}{e_1 e_2 }\left(    v \;\frac{\partial e_2 }{\partial i}
                                          -u  \;\frac{\partial e_1 }{\partial j}  \right) \\
+                                         -u  \;\frac{\partial e_1 }{\partial j}  \right) . \\
   \end{array}
+  }
 …
        -  e_2 u \;\frac{\partial e_3 }{\partial i}
        -  e_1 v \;\frac{\partial e_3 }{\partial j}   \right)
        -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
+       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
     && - \frac{v}{e_1 e_2 }\left(   v  \;\frac{\partial e_2 }{\partial i}
        -u  \;\frac{\partial e_1 }{\partial j}   \right) \\ \\
 …
     && - \,u \left[  \frac{1}{e_1 e_2 e_3} \left(   \frac{\partial(e_2 e_3 \, u)}{\partial i}
        + \frac{\partial(e_1 e_3 \, v)}{\partial j}  \right)
        -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right]
+       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right]
        - \frac{v}{e_1 e_2 }\left(   v   \;\frac{\partial e_2 }{\partial i}
        -u   \;\frac{\partial e_1 }{\partial j}  \right)                  \\
+       -u   \;\frac{\partial e_1 }{\partial j}  \right)     .             \\
+     %
     \intertext {Introducing a more compact form for the divergence of the momentum fluxes,
 …
+  %
     &= \left. {\frac{\partial u }{\partial t}} \right|_s
     &+ \left.  \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)    \right|_s
+    &+ \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s
       + \,u \frac{1}{e_3 } \frac{\partial e_3}{\partial t}
       - \frac{v}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
 …
   \label{apdx:A_sco_Dt_flux}
   \left. \frac{D u}{D t} \right|_s   = \frac{1}{e_3}  \left. \frac{\partial ( e_3\,u)}{\partial t} \right|_s
   + \left.  \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)    \right|_s
+  + \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s
   - \frac{v}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
     -u  \;\frac{\partial e_1 }{\partial j}            \right)
+    -u  \;\frac{\partial e_1 }{\partial j}            \right).
 \end{flalign}
 which is the total time derivative expressed in the curvilinear $s-$coordinate system.
 …
     & =-\frac{1}{\rho_o e_1 }\left[ {\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial p}{\partial s}} \right] \\
     & =-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s +\frac{\sigma_1 }{\rho_o \,e_3 }\left( {-g\;\rho \;e_3 } \right) \\
     &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1
+    &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 .
   \end{split}
 \]
 Applying similar manipulation to the second component and
 replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it comes:
+replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes:
 \begin{equation}
   \label{apdx:A_grad_p_1}
 …
     -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z
     &=-\frac{1}{\rho_o \,e_2 } \left(    \left.               {\frac{\partial p}{\partial j}} \right|_s
       + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) \\
+      + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) . \\
   \end{split}
 \end{equation}
 …
 As in $z$-coordinate,
 the horizontal pressure gradient can be split in two parts following \citet{Marsaleix_al_OM08}.
+the horizontal pressure gradient can be split in two parts following \citet{marsaleix.auclair.ea_OM08}.
 Let defined a density anomaly, $d$, by $d=(\rho - \rho_o)/ \rho_o$,
 and a hydrostatic pressure anomaly, $p_h'$, by $p_h'= g \; \int_z^\eta d \; e_3 \; dk$.
 …
 \[
   \begin{split}
     p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \left(  \rho_o \, d + 1 \right) \; e_3 \; dk   \\
     &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + g \, \int_z^\eta e_3 \; dk
+    p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \rho_o \left( d + 1 \right) \; e_3 \; dk   \\
+    &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + \rho_o g \, \int_z^\eta e_3 \; dk .
   \end{split}
 \]
 …
 \begin{equation}
   \label{apdx:A_pressure}
   p = \rho_o \; p_h' + g \, ( z + \eta )
+  p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z )
 \end{equation}
 and the hydrostatic pressure balance expressed in terms of $p_h'$ and $d$ is:
 \[
   \frac{\partial p_h'}{\partial k} = - d \, g \, e_3
+  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 .
 \]
 Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and
 using the definition of the density anomaly it comes the expression in two parts:
+using the definition of the density anomaly it becomes an expression in two parts:
 \begin{equation}
   \label{apdx:A_grad_p_2}
 …
     -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
     &=-\frac{1}{e_1 } \left(     \left.              {\frac{\partial p_h'}{\partial i}} \right|_s
       + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} \\
+      + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} ,  \\
+             %
     -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z
     &=-\frac{1}{e_2 } \left(    \left.               {\frac{\partial p_h'}{\partial j}} \right|_s
       + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j}\\
+      + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} . \\
   \end{split}
 \end{equation}
 …
     -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right)
     -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i}
     +   D_u^{\vect{U}}  +   F_u^{\vect{U}}
+    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} ,
   \end{multline}
   \begin{multline}
 …
     -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right)
     -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j}
     +  D_v^{\vect{U}}  +   F_v^{\vect{U}}
+    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .
   \end{multline}
 \end{subequations}
 …
     \label{apdx:A_PE_dyn_flux_u}
     \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} =
     \nabla \cdot \left(   {{\rm {\bf U}}\,u}   \right)
+    - \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)
     +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
           -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,v     \\
     -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right)
     -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i}
     +   D_u^{\vect{U}}  +   F_u^{\vect{U}}
+    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} ,
   \end{multline}
   \begin{multline}
     \label{apdx:A_dyn_flux_v}
     \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}=
     -  \nabla \cdot \left(   {{\rm {\bf U}}\,v}   \right)
     +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
+    -  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,v}   \right)
+    -   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
           -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,u     \\
     -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right)
     -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j}
     +  D_v^{\vect{U}}  +   F_v^{\vect{U}}
+    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .
   \end{multline}
 \end{subequations}
 …
 \begin{equation}
   \label{apdx:A_dyn_zph}
   \frac{\partial p_h'}{\partial k} = - d \, g \, e_3
+  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 .
 \end{equation}
 …
   \left[           \frac{\partial }{\partial i} \left( {e_2 \,e_3 \;Tu} \right)
     +   \frac{\partial }{\partial j} \left( {e_1 \,e_3 \;Tv} \right)               \right]       \\
   +  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right)
+  -  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right)
   +  D^{T} +F^{T}
 \end{multline}
 The expression for the advection term is a straight consequence of (A.4),
+The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}),
 the expression of the 3D divergence in the $s-$coordinates established above.

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_B.tex

-                      r10442
+                      r11422
+  {
   \begin{array}{*{20}l}
+    D^T=& \frac{1}{e_1\,e_2\,e_3 }\;\left[ {\ \ \ \ e_2\,e_3\,A^{lT} \;\left.
+          {\frac{\partial }{\partial i}\left( {\frac{1}{e_1}\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_1 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.  \\
+        &\qquad \quad \ \ \ +e_1\,e_3\,A^{lT} \;\left. {\frac{\partial }{\partial j}\left( {\frac{1}{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s -\frac{\sigma_2 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s \\
+        &\qquad \quad \ \ \ + e_1\,e_2\,A^{lT} \;\frac{\partial }{\partial s}\left( {-\frac{\sigma_1 }{e_1 }\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_2 }{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s } \right.
+          \left. {\left. {+\left( {\varepsilon +\sigma_1^2+\sigma_2 ^2} \right)\;\frac{1}{e_3 }\;\frac{\partial T}{\partial s}} \right)\;\;} \right]
+    D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i}  \left. \left[  e_2\,e_3 \, A^{lT}
+                               \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s
+                                       -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
+        &  \quad \  +   \            \left.   \frac{\partial }{\partial j}  \left. \left[  e_1\,e_3 \, A^{lT}
+                               \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s
+                                       -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
+        &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left(
+                     -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s
+                     -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s
+                          +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \;  \right\} .
   \end{array}
+          }
 \end{align*}
 Equation \autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption.
+\autoref{apdx:B2} is obtained from \autoref{apdx:B1} without any additional assumption.
 Indeed, for the special case $k=z$ and thus $e_3 =1$,
 we introduce an arbitrary vertical coordinate $s = s (i,j,z)$ as in \autoref{apdx:A} and
 …
+  {
   \begin{array}{*{20}l}
     \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, it becomes:}
+    \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, this becomes:}
+    %
     & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
+    D^T & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
     & \qquad \qquad \quad -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s -e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\
     & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial z}} \right)\;\;\;} \right] }\\
+    & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }\\
     \\
     &=\frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. {-\frac{\partial }{\partial i}\left( {e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
     & \qquad \qquad \quad \left. {+\frac{e_2 \,\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s} \;\frac{\partial e_3 }{\partial i}}  \right|_s -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s \\
     & \qquad \qquad \quad-e_2 \,\sigma_1 \frac{\partial}{\partial s}\left( {A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right) \\
     & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]}
+    & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} .
   \end{array}
       } \\
 …
   \begin{array}{*{20}l}
+    %
     \intertext{using the same remark as just above, it becomes:}
+    \intertext{Using the same remark as just above, $D^T$ becomes:}
+    %
     &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\
+   D^T &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\
     & \qquad \qquad \quad+\frac{e_1 \,e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}\;\frac{\partial \sigma_1 }{\partial s} - \frac {\sigma_1 }{e_3} A^{lT} \;\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s}\;\frac{\partial T}{\partial s} \\
     & \qquad \qquad \quad-e_2 \left( {A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s +\frac{\partial }{\partial s}\left( {\sigma_1 A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)-\frac{\partial \sigma_1 }{\partial s}\;A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\
     & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }
+    & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] . }
   \end{array}
       } \\
 …
+    %
     \intertext{Since the horizontal scale factors do not depend on the vertical coordinate,
     the last term of the first line and the first term of the last line cancel, while
     the second line reduces to a single vertical derivative, so it becomes:}
+    the two terms on the second line cancel, while
+    the third line reduces to a single vertical derivative, so it becomes:}
+  %
     & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
+    D^T & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
     & \qquad \qquad \quad \shoveright{ \left. {+\frac{\partial }{\partial s}\left( {-e_2 \,\sigma_1 \,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s +A^{lT}\frac{e_1 \,e_2 }{e_3 }\;\left( {\varepsilon +\sigma_1 ^2} \right)\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} \\
+    %
     \intertext{in other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:}
+    \intertext{In other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:}
   \end{array}
   } \\
 …
 the ($i$,$j$,$k$) curvilinear coordinate system in which
 the equations of the ocean circulation model are formulated,
 takes the following form \citep{Redi_JPO82}:
+takes the following form \citep{redi_JPO82}:
 \begin{equation}
 …
   \left[ {{
         \begin{array}{*{20}c}
           {1+a_1 ^2} \hfill & {-a_1 a_2 } \hfill & {-a_1 } \hfill \\
           {-a_1 a_2 } \hfill & {1+a_2 ^2} \hfill & {-a_2 } \hfill \\
           {-a_1 } \hfill & {-a_2 } \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\
+          {1+a_2 ^2 +\varepsilon a_1 ^2} \hfill & {-a_1 a_2 (1-\varepsilon)} \hfill & {-a_1 (1-\varepsilon) } \hfill \\
+          {-a_1 a_2 (1-\varepsilon) } \hfill & {1+a_1 ^2 +\varepsilon a_2 ^2} \hfill & {-a_2 (1-\varepsilon)} \hfill \\
+          {-a_1 (1-\varepsilon)} \hfill & {-a_2 (1-\varepsilon)} \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\
         \end{array}
       }} \right]
 \end{equation}
+where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials:
+where ($a_1$, $a_2$) are $(-1) \times$ the isopycnal slopes in
+($\textbf{i}$, $\textbf{j}$) directions, relative to geopotentials (or
+equivalently the slopes of the geopotential surfaces in the isopycnal
+coordinate framework):
 \[
   a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
 …
   \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
 \]
+In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean,
+so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}:
+and, as before, $\epsilon = A^{vT} / A^{lT}$.
+In practice, $\epsilon$ is small and isopycnal slopes are generally less than $10^{-2}$ in the ocean,
+so $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{cox_OM87}. Keeping leading order terms\footnote{Apart from the (1,0)
+and (0,1) elements which are set to zero. See \citet{griffies_bk04}, section 14.1.4.1 for a discussion of this point.}:
 \begin{subequations}
   \label{apdx:B4}
 …
 The isopycnal diffusion operator \autoref{apdx:B4},
 \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square.
+The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes.
 Let us demonstrate the second one:
+As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero
+(as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one:
 \[
   \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
 …
+  {
   \begin{array}{*{20}l}
     \nabla T\;.\left( {{\rm {\bf A}}_{\rm {\bf I}} \nabla T}
+    \nabla T\;.\left( {{\mathrm {\mathbf A}}_{\mathrm {\mathbf I}} \nabla T}
     \right)&=A^{lT}\left[ {\left( {\frac{\partial T}{\partial i}} \right)^2-2a_1
              \frac{\partial T}{\partial i}\frac{\partial T}{\partial k}+\left(
 …
              j}-a_2 \frac{\partial T}{\partial k}} \right)^2}
              +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\
            & \geq 0
+           & \geq 0 .
   \end{array}
+             }
 …
 \end{equation}
+where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions,
+relative to $s$-coordinate surfaces:
+where ($r_1$, $r_2$) are $(-1)\times$ the isopycnal slopes in ($\textbf{i}$, $\textbf{j}$) directions,
+relative to $s$-coordinate surfaces (or equivalently the slopes of the
+$s$-coordinate surfaces in the isopycnal coordinate framework):
 \[
   r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1}
 …
 \]
 To prove \autoref{apdx:B5} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious.
+To prove \autoref{apdx:B_ldfiso_s} by direct re-expression of \autoref{apdx:B_ldfiso} is straightforward, but laborious.
 An easier way is first to note (by reversing the derivation of \autoref{sec:B_2} from \autoref{sec:B_1} ) that
 the weak-slope operator may be \emph{exactly} reexpressed in non-orthogonal $i,j,\rho$-coordinates as
 …
 the (rotated,orthogonal) isoneutral axes to the non-orthogonal $i,j,\rho$-coordinates.
 The further transformation into $i,j,s$-coordinates is exact, whatever the steepness of the $s$-surfaces,
 in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates,
+in the same way as the transformation of horizontal/vertical Laplacian diffusion in $z$-coordinates in
 \autoref{sec:B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
 …
 \label{sec:B_3}
 The second order momentum diffusion operator (Laplacian) in the $z$-coordinate is found by
+The second order momentum diffusion operator (Laplacian) in $z$-coordinates is found by
 applying \autoref{eq:PE_lap_vector}, the expression for the Laplacian of a vector,
 to the horizontal velocity vector:
 …
 \end{align*}
 Using \autoref{eq:PE_div}, the definition of the horizontal divergence,
 the third componant of the second vector is obviously zero and thus :
+the third component of the second vector is obviously zero and thus :
 \[
   \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right)
+  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) .
 \]
 …
   - \nabla _h \times \left( {A^{lm}\;\zeta \;{\textbf{k}}} \right)
   + \frac{1}{e_3 }\frac{\partial }{\partial k}\left( {\frac{A^{vm}\;}{e_3 }
       \frac{\partial {\rm {\bf U}}_h }{\partial k}} \right) \\
+      \frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}} \right) , \\
 \end{equation}
 that is, in expanded form:
 …
   & = \frac{1}{e_1} \frac{\partial \left( {A^{lm}\chi   } \right)}{\partial i}
     -\frac{1}{e_2} \frac{\partial \left( {A^{lm}\zeta } \right)}{\partial j}
     +\frac{1}{e_3} \frac{\partial u}{\partial k}      \\
+    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial u}{\partial k} \right)   ,   \\
   D^{\textbf{U}}_v
   & = \frac{1}{e_2 }\frac{\partial \left( {A^{lm}\chi   } \right)}{\partial j}
     +\frac{1}{e_1 }\frac{\partial \left( {A^{lm}\zeta } \right)}{\partial i}
     +\frac{1}{e_3} \frac{\partial v}{\partial k}
+    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial v}{\partial k} \right) .
 \end{align*}

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_E.tex

-                      r10442
+                      r11422
 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error
 \citep{Shchepetkin_McWilliams_OM05}.
 The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 It is not a \emph{positive} scheme meaning false extrema are permitted but
 …
 the second term which is the diffusive part of the scheme, is evaluated using the \textit{before} velocity
 (forward in time).
 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
+This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
 UBS and QUICK schemes only differ by one coefficient.
 Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Substituting 1/6 with 1/8 in (\autoref{eq:tra_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 …
 $\tau_w^{ubs}$ will be evaluated using either \textit{(a)} a centered $2^{nd}$ order scheme,
 or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative parabolic splines following
 \citet{Shchepetkin_McWilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
+\citet{shchepetkin.mcwilliams_OM05} implementation of UBS in ROMS, or \textit{(d)} an UBS.
 The $3^{rd}$ case has dispersion properties similar to an eight-order accurate conventional scheme.
 …
 \subsection{Griffies iso-neutral diffusion operator}
 Let try to define a scheme that get its inspiration from the \citet{Griffies_al_JPO98} scheme,
+Let try to define a scheme that get its inspiration from the \citet{griffies.gnanadesikan.ea_JPO98} scheme,
 but is formulated within the \NEMO framework
 (\ie using scale factors rather than grid-size and having a position of $T$-points that
 …
 Nevertheless, this technique works fine for $T$ and $S$ as they are active tracers
 (\ie they enter the computation of density), but it does not work for a passive tracer.
 \citep{Griffies_al_JPO98} introduce a different way to discretise the off-diagonal terms that
+\citep{griffies.gnanadesikan.ea_JPO98} introduce a different way to discretise the off-diagonal terms that
 nicely solve the problem.
 The idea is to get rid of combinations of an averaged in one direction combined with
 …
 \begin{figure}[!ht]
   \begin{center}
     \includegraphics[width=0.70\textwidth]{Fig_ISO_triad}
+    \includegraphics[width=\textwidth]{Fig_ISO_triad}
     \caption{
       \protect\label{fig:ISO_triad}
 …
 \]
 \citep{Griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
+\citep{griffies_JPO98} introduces another way to implement the eddy induced advection, the so-called skew form.
 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
 For example in the (\textbf{i},\textbf{k}) plane, the tracer advective fluxes can be transformed as follows:
 …
 The horizontal component reduces to the one use for an horizontal laplacian operator and
 the vertical one keeps the same complexity, but not more.
 This property has been used to reduce the computational time \citep{Griffies_JPO98},
+This property has been used to reduce the computational time \citep{griffies_JPO98},
 but it is not of practical use as usually $A \neq A_e$.
 Nevertheless this property can be used to choose a discret form of \autoref{eq:eiv_skew_continuous} which

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_iso.tex

-                      r10442
+                      r11422
 \newcommand{\rML}[1][i]{\ensuremath{_{\mathrm{ML}\,#1}}}
 \newcommand{\rMLt}[1][i]{\tilde{r}_{\mathrm{ML}\,#1}}
+\newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}}}
+%% Move to ../../global/new_cmds.tex to avoid error with \listoffigures
+%\newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}}
 \newcommand{\triadd}[5]{\ensuremath{{}_{#1}^{#2}{\mathbb{#3}}_{#4}^{\,#5}}}
 \newcommand{\triadt}[5]{\ensuremath{{}_{#1}^{#2}{\tilde{\mathbb{#3}}}_{#4}^{\,#5}}}
 …
   the vertical skew flux is further reduced to ensure no vertical buoyancy flux,
   giving an almost pure horizontal diffusive tracer flux within the mixed layer.
   This is similar to the tapering suggested by \citet{Gerdes1991}. See \autoref{subsec:Gerdes-taper}
+  This is similar to the tapering suggested by \citet{gerdes.koberle.ea_CD91}. See \autoref{subsec:Gerdes-taper}
 \item[\np{ln\_botmix\_triad}]
   See \autoref{sec:iso_bdry}.
 …
 \label{sec:iso}
 We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98},
+We have implemented into \NEMO a scheme inspired by \citet{griffies.gnanadesikan.ea_JPO98},
 but formulated within the \NEMO framework, using scale factors rather than grid-sizes.
 …
 \subsection{Expression of the skew-flux in terms of triad slopes}
 \citep{Griffies_al_JPO98} introduce a different discretization of the off-diagonal terms that
+\citep{griffies.gnanadesikan.ea_JPO98} introduce a different discretization of the off-diagonal terms that
 nicely solves the problem.
 % Instead of multiplying the mean slope calculated at the $u$-point by
 …
 \begin{figure}[tb]
   \begin{center}
     \includegraphics[width=1.05\textwidth]{Fig_GRIFF_triad_fluxes}
+    \includegraphics[width=\textwidth]{Fig_GRIFF_triad_fluxes}
     \caption{
       \protect\label{fig:ISO_triad}
 …
 \begin{figure}[tb]
   \begin{center}
     \includegraphics[width=0.80\textwidth]{Fig_GRIFF_qcells}
+    \includegraphics[width=\textwidth]{Fig_GRIFF_qcells}
     \caption{
       \protect\label{fig:qcells}
 …
 To complete the discretization we now need only specify the triad volumes $_i^k\mathbb{V}_{i_p}^{k_p}$.
 \citet{Griffies_al_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
+\citet{griffies.gnanadesikan.ea_JPO98} identifies these $_i^k\mathbb{V}_{i_p}^{k_p}$ as the volumes of the quarter cells,
 defined in terms of the distances between $T$, $u$,$f$ and $w$-points.
 This is the natural discretization of \autoref{eq:cts-var}.
 …
 \begin{figure}[h]
   \begin{center}
     \includegraphics[width=0.60\textwidth]{Fig_GRIFF_bdry_triads}
+    \includegraphics[width=\textwidth]{Fig_GRIFF_bdry_triads}
     \caption{
       \protect\label{fig:bdry_triads}
 …
 As discussed in \autoref{subsec:LDF_slp_iso},
 iso-neutral slopes relative to geopotentials must be bounded everywhere,
 both for consistency with the small-slope approximation and for numerical stability \citep{Cox1987, Griffies_Bk04}.
+both for consistency with the small-slope approximation and for numerical stability \citep{cox_OM87, griffies_bk04}.
 The bound chosen in \NEMO is applied to each component of the slope separately and
 has a value of $1/100$ in the ocean interior.
 …
 \[
   % \label{eq:iso_tensor_ML}
   D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
+  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
   \mbox{with}\quad \;\;\Re =\left( {{
         \begin{array}{*{20}c}
 …
     (\eg the green triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}
   % }
   \includegraphics[width=0.60\textwidth]{Fig_GRIFF_MLB_triads}
+  \includegraphics[width=\textwidth]{Fig_GRIFF_MLB_triads}
 \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 \[
   % \label{eq:iso_tensor_ML2}
   D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
+  D^{lT}=\nabla {\mathrm {\mathbf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
   \mbox{with}\quad \;\;\Re =\left( {{
         \begin{array}{*{20}c}
 …
 \footnote{
   To ensure good behaviour where horizontal density gradients are weak,
   we in fact follow \citet{Gerdes1991} and
+  we in fact follow \citet{gerdes.koberle.ea_CD91} and
   set $\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.
+}
 …
 This approach is similar to multiplying the iso-neutral diffusion coefficient by
 $\tilde{r}_{\mathrm{max}}^{-2}\tilde{r}_i^{-2}$ for steep slopes,
 as suggested by \citet{Gerdes1991} (see also \citet{Griffies_Bk04}).
+as suggested by \citet{gerdes.koberle.ea_CD91} (see also \citet{griffies_bk04}).
 Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
 …
 However, when \np{ln\_traldf\_triad} is set true,
 \NEMO instead implements eddy induced advection according to the so-called skew form \citep{Griffies_JPO98}.
+\NEMO instead implements eddy induced advection according to the so-called skew form \citep{griffies_JPO98}.
 It is based on a transformation of the advective fluxes using the non-divergent nature of the eddy induced velocity.
 For example in the (\textbf{i},\textbf{k}) plane,
 …
 it is equivalent to a horizontal eiv (eddy-induced velocity) that is uniform within the mixed layer
 \autoref{eq:eiv_v}.
 This ensures that the eiv velocities do not restratify the mixed layer \citep{Treguier1997,Danabasoglu_al_2008}.
+This ensures that the eiv velocities do not restratify the mixed layer \citep{treguier.held.ea_JPO97,danabasoglu.ferrari.ea_JC08}.
 Equivantly, in terms of the skew-flux formulation we use here,
 the linear slope tapering within the mixed-layer gives a linearly varying vertical flux,
 …
 $uw$ (integer +1/2 $i$, integer $j$, integer +1/2 $k$) and $vw$ (integer $i$, integer +1/2 $j$, integer +1/2 $k$)
 points (see Table \autoref{tab:cell}) respectively.
 We follow \citep{Griffies_Bk04} and calculate the streamfunction at a given $uw$-point from
+We follow \citep{griffies_bk04} and calculate the streamfunction at a given $uw$-point from
 the surrounding four triads according to:
 \[

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_ASM.tex

-                      r10442
+                      r11422
 \label{chap:ASM}
+Authors: D. Lea,  M. Martin, K. Mogensen, A. Weaver, ...   % do we keep
+\minitoc
+\minitoc
+\vfill
+\begin{figure}[b]
+\subsubsection*{Changes record}
+\begin{tabular}{l||l|m{0.65\linewidth}}
+    Release   & Author        & Modifications \\
+    {\em 4.0} & {\em D. J. Lea} & {\em \NEMO 4.0 updates}  \\
+    {\em 3.4} & {\em D. J. Lea, M. Martin, K. Mogensen, A. Weaver} & {\em Initial version}  \\
+\end{tabular}
+\end{figure}
 \newpage
 The ASM code adds the functionality to apply increments to the model variables: temperature, salinity,
 sea surface height, velocity and sea ice concentration.
+sea surface height, velocity and sea ice concentration.
 These are read into the model from a NetCDF file which may be produced by separate data assimilation code.
 The code can also output model background fields which are used as an input to data assimilation code.
 …
 it may be preferable to introduce the increment gradually into the ocean model in order to
 minimize spurious adjustment processes.
 This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}.
+This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}.
 IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI.
 IAU is used when \np{ln\_asmiau} is set to true.
 With IAU, the model state trajectory ${\bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$)
+With IAU, the model state trajectory ${\mathbf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$)
 is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as
 additional tendency terms to the prognostic equations:
 \begin{align*}
   % \label{eq:wa_traj_iau}
   {\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\bf x}^{a}
+  {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a}
 \end{align*}
 where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\bf x}^{a}$ defined such that
+where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\mathbf x}^{a}$ defined such that
 $\sum_{i=1}^{N} F_{i}=1$.
 ${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied.
+${\mathbf x}^b$ denotes the model initial state and ${\mathbf x}^a$ is the model state after the increments are applied.
 To control the adjustment time of the model to the increment,
 the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$,
 …
 Typically the increments are spread evenly over the full window.
 In addition, two different weighting functions have been implemented.
 The first function employs constant weights,
+The first function (namelist option \np{niaufn} = 0) employs constant weights,
 \begin{align}
   \label{eq:F1_i}
 …
   =\left\{
   \begin{array}{ll}
      &    {\rm if} \; \; \; t_{i} < t_{m}                \\
 /M &    {\rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\
      &    {\rm if} \; \; \; t_{i} > t_{n}
+     &    {\mathrm if} \; \; \; t_{i} < t_{m}                \\
+/M &    {\mathrm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\
+     &    {\mathrm if} \; \; \; t_{i} > t_{n}
   \end{array}
             \right.
+            \right.
 \end{align}
 where $M = m-n$.
 The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window,
+The second function (namelist option \np{niaufn} = 1) employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window,
 with the weighting reduced linearly to a small value at the window end-points:
 \begin{align}
 …
   =\left\{
   \begin{array}{ll}
                            &    {\rm if} \; \; \; t_{i}       <     t_{m}                        \\
     \alpha \, i               &    {\rm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\
     \alpha \, (M - i +1) &    {\rm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\
                             &   {\rm if} \; \; \; t_{i}        >    t_{n}
+                           &    {\mathrm if} \; \; \; t_{i}       <     t_{m}                        \\
+    \alpha \, i               &    {\mathrm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\
+    \alpha \, (M - i +1) &    {\mathrm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\
+                            &   {\mathrm if} \; \; \; t_{i}        >    t_{n}
   \end{array}
                                    \right.
 \end{align}
 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.
+where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.
 The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from
 one assimilation cycle to the next than that described by \autoref{eq:F1_i}.
 …
 \label{sec:ASM_div_dmp}
+The velocity increments may be initialized by the iterative application of a divergence damping operator.
+In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by:
+It is quite challenging for data assimilation systems to provide non-divergent velocity increments.
+Applying divergent velocity increments will likely cause spurious vertical velocities in the model. This section describes a method to take velocity increments provided to \NEMO ($u^0_I$ and $v^0_I$) and adjust them by the iterative application of a divergence damping operator. The method is also described in \citet{dobricic.pinardi.ea_OS07}.
+In iteration step $n$ (starting at $n=1$) new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by:
 \begin{equation}
   \label{eq:asm_dmp}
 …
   \right.,
 \end{equation}
+where
+where the divergence is defined as
 \[
   % \label{eq:asm_div}
 …
       +\delta_j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right).
 \]
+By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration,
+By the application of \autoref{eq:asm_dmp} the divergence is filtered in each iteration,
 and the vorticity is left unchanged.
 In the presence of coastal boundaries with zero velocity increments perpendicular to the coast
 …
 This type of the initialisation reduces the vertical velocity magnitude and
 alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration
 \citep{Talagrand_JAS72, Dobricic_al_OS07}.
+\citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}.
 Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$.
 The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist
 a value greater than zero.
+By choosing this value to be of the order of 100 the increments in
+the vertical velocity will be significantly reduced.
+This specifies the number of iterations of the divergence damping. Setting a value of the order of 100 will result in a significant reduction in the vertical velocity induced by the increments.
 …
 \label{sec:ASM_details}
 Here we show an example \ngn{namasm} namelist and the header of an example assimilation increments file on
+Here we show an example \ngn{nam\_asminc} namelist and the header of an example assimilation increments file on
 the ORCA2 grid.
 %------------------------------------------namasm-----------------------------------------------------
+%------------------------------------------nam_asminc-----------------------------------------------------
+%
 \nlst{nam_asminc}

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_CONFIG.tex

-                      r10442
+                      r11422
 \label{sec:CFG_intro}
 The purpose of this part of the manual is to introduce the \NEMO reference configurations.
+The purpose of this part of the manual is to introduce the NEMO reference configurations.
 These configurations are offered as means to explore various numerical and physical options,
 thus allowing the user to verify that the code is performing in a manner consistent with that we are running.
 …
 The reference configurations also provide a sense for some of the options available in the code,
 though by no means are all options exercised in the reference configurations.
+Configuration is defined manually through the \textit{namcfg} namelist variables.
 %------------------------------------------namcfg----------------------------------------------------
 …
 % 1D model configuration
 % ================================================================
+\section{C1D: 1D Water column model (\protect\key{c1d}) }
+\section[C1D: 1D Water column model (\texttt{\textbf{key\_c1d}})]
+{C1D: 1D Water column model (\protect\key{c1d})}
 \label{sec:CFG_c1d}
+BE careful: to be re-written according to suppression of jpizoom and jpjzoom !!!!
+The 1D model option simulates a stand alone water column within the 3D \NEMO system.
+The 1D model option simulates a stand alone water column within the 3D NEMO system.
 It can be applied to the ocean alone or to the ocean-ice system and can include passive tracers or a biogeochemical model.
 It is set up by defining the position of the 1D water column in the grid
 (see \textit{CONFIG/SHARED/namelist\_ref} ).
+(see \textit{cfgs/SHARED/namelist\_ref}).
 The 1D model is a very useful tool
 \textit{(a)} to learn about the physics and numerical treatment of vertical mixing processes;
 …
 \textit{(d)} to produce extra diagnostics, without the large memory requirement of the full 3D model.
+The methodology is based on the use of the zoom functionality over the smallest possible domain:
+a 3x3 domain centered on the grid point of interest, with some extra routines.
+There is no need to define a new mesh, bathymetry, initial state or forcing,
+since the 1D model will use those of the configuration it is a zoom of.
+The chosen grid point is set in \textit{\ngn{namcfg}} namelist by
+setting the \np{jpizoom} and \np{jpjzoom} parameters to the indices of the location of the chosen grid point.
+The methodology is based on the configuration of the smallest possible domain:
+a 3x3 domain with 75 vertical levels.
 The 1D model has some specifies. First, all the horizontal derivatives are assumed to be zero,
 and second, the two components of the velocity are moved on a $T$-point.
 Therefore, defining \key{c1d} changes five main things in the code behaviour:
+Therefore, defining \key{c1d} changes some things in the code behaviour:
 \begin{description}
 \item[(1)]
-  the lateral boundary condition routine (\rou{lbc\_lnk}) set the value of the central column of
-  the 3x3 domain is imposed over the whole domain;
-\item[(3)]
-  a call to \rou{lbc\_lnk} is systematically done when reading input data (\ie in \mdl{iom});
-\item[(3)]
   a simplified \rou{stp} routine is used (\rou{stp\_c1d}, see \mdl{step\_c1d} module) in which
   both lateral tendancy terms and lateral physics are not called;
 \item[(4)]
+\item[(2)]
   the vertical velocity is zero
   (so far, no attempt at introducing a Ekman pumping velocity has been made);
 \item[(5)]
+\item[(3)]
   a simplified treatment of the Coriolis term is performed as $U$- and $V$-points are the same
   (see \mdl{dyncor\_c1d}).
 \end{description}
 All the relevant \textit{\_c1d} modules can be found in the NEMOGCM/NEMO/OPA\_SRC/C1D directory of
 the \NEMO distribution.
+All the relevant \textit{\_c1d} modules can be found in the src/OCE/C1D directory of
+the NEMO distribution.
 % to be added:  a test case on the yearlong Ocean Weather Station (OWS) Papa dataset of Martin (1985)
 …
 The ORCA family is a series of global ocean configurations that are run together with
+the LIM sea-ice model (ORCA-LIM) and possibly with PISCES biogeochemical model (ORCA-LIM-PISCES),
+using various resolutions.
+An appropriate namelist is available in \path{CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg} for ORCA2.
+the SI3 model (ORCA-ICE) and possibly with PISCES biogeochemical model (ORCA-ICE-PISCES).
+An appropriate namelist is available in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg} for ORCA2.
 The domain of ORCA2 configuration is defined in \ifile{ORCA\_R2\_zps\_domcfg} file,
+this file is available in tar file in the wiki of NEMO: \\
+https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ReferenceConfigurations/ORCA2\_LIM3\_PISCES \\
+this file is available in tar file on the NEMO community zenodo platform: \\
+https://doi.org/10.5281/zenodo.2640723
 In this namelist\_cfg the name of domain input file is set in \ngn{namcfg} block of namelist.
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=0.98\textwidth]{Fig_ORCA_NH_mesh}
+    \includegraphics[width=\textwidth]{Fig_ORCA_NH_mesh}
     \caption{
       \protect\label{fig:MISC_ORCA_msh}
 …
       The two "north pole" are the foci of a series of embedded ellipses (blue curves) which
       are determined analytically and form the i-lines of the ORCA mesh (pseudo latitudes).
       Then, following \citet{Madec_Imbard_CD96}, the normal to the series of ellipses (red curves) is computed which
+      Then, following \citet{madec.imbard_CD96}, the normal to the series of ellipses (red curves) is computed which
       provides the j-lines of the mesh (pseudo longitudes).
+    }
 …
 \label{subsec:CFG_orca_grid}
 The ORCA grid is a tripolar is based on the semi-analytical method of \citet{Madec_Imbard_CD96}.
+The ORCA grid is a tripolar grid based on the semi-analytical method of \citet{madec.imbard_CD96}.
 It allows to construct a global orthogonal curvilinear ocean mesh which has no singularity point inside
 the computational domain since two north mesh poles are introduced and placed on lands.
 …
 \begin{figure}[!tbp]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_ORCA_NH_msh05_e1_e2}
     \includegraphics[width=0.80\textwidth]{Fig_ORCA_aniso}
+    \includegraphics[width=\textwidth]{Fig_ORCA_NH_msh05_e1_e2}
+    \includegraphics[width=\textwidth]{Fig_ORCA_aniso}
     \caption {
       \protect\label{fig:MISC_ORCA_e1e2}
 …
 % -------------------------------------------------------------------------------------------------------------
 %       ORCA-LIM(-PISCES) configurations
+%       ORCA-ICE(-PISCES) configurations
 % -------------------------------------------------------------------------------------------------------------
 \subsection{ORCA pre-defined resolution}
 …
 The ORCA\_R2 configuration has the following specificity: starting from a 2\deg~ORCA mesh,
 local mesh refinements were applied to the Mediterranean, Red, Black and Caspian Seas,
 so that the resolution is 1\deg \time 1\deg there.
+so that the resolution is 1\deg~ there.
 A local transformation were also applied with in the Tropics in order to refine the meridional resolution up to
 .5\deg at the Equator.
+.5\deg~ at the Equator.
 The ORCA\_R1 configuration has only a local tropical transformation to refine the meridional resolution up to
 …
 For ORCA\_R1 and R025, setting the configuration key to 75 allows to use 75 vertical levels, otherwise 46 are used.
 In the other ORCA configurations, 31 levels are used
 (see \autoref{tab:orca_zgr} %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}).
 Only the ORCA\_R2 is provided with all its input files in the \NEMO distribution.
 It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of
 climate change (Marti et al., 2009).
 It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs)
 documented in \citet{Griffies_al_OM09}.
+(see \autoref{tab:orca_zgr}). %\sfcomment{HERE I need to put new table for ORCA2 values} and \autoref{fig:zgr}).
+Only the ORCA\_R2 is provided with all its input files in the NEMO distribution.
+%It is very similar to that used as part of the climate model developed at IPSL for the 4th IPCC assessment of
+%climate change (Marti et al., 2009).
+%It is also the basis for the \NEMO contribution to the Coordinate Ocean-ice Reference Experiments (COREs)
+%documented in \citet{griffies.biastoch.ea_OM09}.
 This version of ORCA\_R2 has 31 levels in the vertical, with the highest resolution (10m) in the upper 150m
 (see \autoref{tab:orca_zgr} and \autoref{fig:zgr}).
 The bottom topography and the coastlines are derived from the global atlas of Smith and Sandwell (1997).
 The default forcing uses the boundary forcing from \citet{Large_Yeager_Rep04} (see \autoref{subsec:SBC_blk_core}),
+The default forcing uses the boundary forcing from \citet{large.yeager_rpt04} (see \autoref{subsec:SBC_blk_core}),
 which was developed for the purpose of running global coupled ocean-ice simulations without
 an interactive atmosphere.
 This \citet{Large_Yeager_Rep04} dataset is available through
+This \citet{large.yeager_rpt04} dataset is available through
 the \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.
 The "normal year" of \citet{Large_Yeager_Rep04} has been chosen of the \NEMO distribution since release v3.3.
 ORCA\_R2 pre-defined configuration can also be run with an AGRIF zoom over the Agulhas current area
+(\key{agrif} defined) and, by setting the appropriate variables, see \path{CONFIG/SHARED/namelist_ref}.
+The "normal year" of \citet{large.yeager_rpt04} has been chosen of the NEMO distribution since release v3.3.
+ORCA\_R2 pre-defined configuration can also be run with multiply online nested zooms (\ie with AGRIF, \key{agrif} defined). This is available as the AGRIF\_DEMO configuration that can be found in the \path{cfgs/AGRIF_DEMO/} directory.
 A regional Arctic or peri-Antarctic configuration is extracted from an ORCA\_R2 or R05 configurations using
 sponge layers at open boundaries.
 …
 %       GYRE family: double gyre basin
 % -------------------------------------------------------------------------------------------------------------
 \section{GYRE family: double gyre basin }
+\section{GYRE family: double gyre basin}
 \label{sec:CFG_gyre}
 The GYRE configuration \citep{Levy_al_OM10} has been built to
+The GYRE configuration \citep{levy.klein.ea_OM10} has been built to
 simulate the seasonal cycle of a double-gyre box model.
 It consists in an idealized domain similar to that used in the studies of \citet{Drijfhout_JPO94} and
 \citet{Hazeleger_Drijfhout_JPO98, Hazeleger_Drijfhout_JPO99, Hazeleger_Drijfhout_JGR00, Hazeleger_Drijfhout_JPO00},
+It consists in an idealized domain similar to that used in the studies of \citet{drijfhout_JPO94} and
+\citet{hazeleger.drijfhout_JPO98, hazeleger.drijfhout_JPO99, hazeleger.drijfhout_JGR00, hazeleger.drijfhout_JPO00},
 over which an analytical seasonal forcing is applied.
 This allows to investigate the spontaneous generation of a large number of interacting, transient mesoscale eddies
 and their contribution to the large scale circulation.
+The GYRE configuration run together with the PISCES biogeochemical model (GYRE-PISCES).
 The domain geometry is a closed rectangular basin on the $\beta$-plane centred at $\sim$ 30\deg{N} and
 rotated by 45\deg, 3180~km long, 2120~km wide and 4~km deep (\autoref{fig:MISC_strait_hand}).
 …
 The configuration is meant to represent an idealized North Atlantic or North Pacific basin.
 The circulation is forced by analytical profiles of wind and buoyancy fluxes.
 The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{Levy_al_OM10}.
+The applied forcings vary seasonally in a sinusoidal manner between winter and summer extrema \citep{levy.klein.ea_OM10}.
 The wind stress is zonal and its curl changes sign at 22\deg{N} and 36\deg{N}.
 It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain and
 …
 The GYRE configuration is set like an analytical configuration.
 Through \np{ln\_read\_cfg}\forcode{ = .false.} in \textit{namcfg} namelist defined in
 the reference configuration \path{CONFIG/GYRE/EXP00/namelist_cfg}
+the reference configuration \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg}
 analytical definition of grid in GYRE is done in usrdef\_hrg, usrdef\_zgr routines.
 Its horizontal resolution (and thus the size of the domain) is determined by
 setting \np{nn\_GYRE} in \ngn{namusr\_def}: \\
 \np{jpiglo} $= 30 \times$ \np{nn\_GYRE} + 2   \\
 \np{jpjglo} $= 20 \times$ \np{nn\_GYRE} + 2   \\
 Obviously, the namelist parameters have to be adjusted to the chosen resolution,
 see the Configurations pages on the NEMO web site (Using NEMO\/Configurations).
+see the Configurations pages on the NEMO web site (NEMO Configurations).
 In the vertical, GYRE uses the default 30 ocean levels (\jp{jpk}\forcode{ = 31}) (\autoref{fig:zgr}).
 …
 even though the physical integrity of the solution can be compromised.
 Benchmark is activate via \np{ln\_bench}\forcode{ = .true.} in \ngn{namusr\_def} in
 namelist \path{CONFIG/GYRE/EXP00/namelist_cfg}.
+namelist \path{cfgs/GYRE_PISCES/EXPREF/namelist_cfg}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_GYRE}
+    \includegraphics[width=\textwidth]{Fig_GYRE}
     \caption{
       \protect\label{fig:GYRE}
       Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54.
       From \citet{Levy_al_OM10}.
+      From \citet{levy.klein.ea_OM10}.
+    }
   \end{center}
 …
 The AMM, Atlantic Margins Model, is a regional model covering the Northwest European Shelf domain on
 a regular lat-lon grid at approximately 12km horizontal resolution.
 The appropriate \textit{\&namcfg} namelist  is available in \textit{CONFIG/AMM12/EXP00/namelist\_cfg}.
+The appropriate \textit{\&namcfg} namelist  is available in \textit{cfgs/AMM12/EXPREF/namelist\_cfg}.
 It is used to build the correct dimensions of the AMM domain.
 This configuration tests several features of NEMO functionality specific to the shelf seas.
+In particular, the AMM uses $S$-coordinates in the vertical rather than $z$-coordinates and
+is forced with tidal lateral boundary conditions using a flather boundary condition from the BDY module.
+The AMM configuration uses the GLS (\key{zdfgls}) turbulence scheme,
+the VVL non-linear free surface(\key{vvl}) and time-splitting (\key{dynspg\_ts}).
+In particular, the AMM uses $s$-coordinates in the vertical rather than $z$-coordinates and
+is forced with tidal lateral boundary conditions using a Flather boundary condition from the BDY module.
+Also specific to the AMM configuration is the use of the GLS turbulence scheme (\np{ln\_zdfgls} \forcode{= .true.}).
 In addition to the tidal boundary condition the model may also take open boundary conditions from

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIA.tex

-                      r10509
+                      r11422
 \minitoc
+\vfill
+\begin{figure}[b]
+\subsubsection*{Changes record}
+\begin{tabular}{l||l|m{0.65\linewidth}}
+    Release   & Author        & Modifications \\
+    {\em 4.0} & {\em Mirek Andrejczuk, Massimiliano Drudi} & {\em }  \\
+    {\em }      & {\em Dorotea Iovino, Nicolas Martin} & {\em }  \\
+    {\em 3.6} & {\em Gurvan Madec, Sebastien Masson } & {\em }  \\
+    {\em 3.4} & {\em Gurvan Madec, Rachid Benshila, Andrew Coward } & {\em }  \\
+    {\em }      & {\em Christian Ethe, Sebastien Masson } & {\em }  \\
+\end{tabular}
+\end{figure}
 \newpage
 …
 %       Old Model Output
 % ================================================================
 \section{Old model output (default)}
+\section{Model output}
 \label{sec:DIA_io_old}
 …
 the same run performed in one step.
 It should be noted that this requires that the restart file contains two consecutive time steps for
+all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that
+is to read it (in particular, 32 bits binary IEEE format must not be used for this file).
+all the prognostic variables.
 The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs.
 …
 the writing work is distributed over the processors in massively parallel computing.
 A complete description of the use of this I/O server is presented in the next section.
-By default, \key{iomput} is not defined,
-NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation.
-However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2,
-many diagnostic options have been added presuming the use of \key{iomput}.
-The usefulness of the default IOIPSL-based option is expected to reduce with each new release.
-If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and
-contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of
-nn\_write time-steps (namelist parameter).
 %\gmcomment{                    % start of gmcomment
 …
 in a very easy way.
 All details of iomput functionalities are listed in the following subsections.
+Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml},
+\path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\
+Examples of the XML files that control the outputs can be found in:
+\path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef.xml},
+\path{cfgs/SHARED/field_def_nemo-oce.xml},
+\path{cfgs/SHARED/field_def_nemo-pisces.xml},
+\path{cfgs/SHARED/field_def_nemo-ice.xml} and \path{cfgs/SHARED/domain_def_nemo.xml}. \\
 The second functionality targets output performance when running in parallel (\key{mpp\_mpi}).
 …
 In version 3.6, the iom\_put interface depends on
 an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0}
+an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-2.5}{XIOS-2.5}
 (use of revision 618 or higher is required).
 This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to
 …
 \xmlline|<variable id="using_server" type="bool"></variable>|
 The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode}
 (as a library) [{\tt> false <}] or in \textit{detached mode}
 (as an external executable on N additional, dedicated cpus) [{\tt > true <}].
+The {\ttfamily using\_server} setting determines whether or not the server will be used in \textit{attached mode}
+(as a library) [{\ttfamily> false <}] or in \textit{detached mode}
+(as an external executable on N additional, dedicated cpus) [{\ttfamily > true <}].
 The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications.
 The type of each file can be either ''multiple\_file'' or ''one\_file''.
 …
 \subsubsection{Control of XIOS: the context in iodef.xml}
 As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml.
+As well as the {\ttfamily using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml.
 See the XML basics section below for more details on XML syntax and rules.
 …
 See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance.
 NEMO will need to link to the compiled XIOS library.
 The \href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS}
 {XIOS with NEMO} guide provides an example illustration of how this can be achieved.
+The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/install.html#extract-and-install-xios}
+{Extract and install XIOS} guide provides an example illustration of how this can be achieved.
 \subsubsection{Add your own outputs}
 …
 \begin{enumerate}
 \item[1.]
   in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to output a 2D or 3D array.
+  in NEMO code, add a \forcode{CALL iom_put( 'identifier', array )} where you want to output a 2D or 3D array.
 \item[2.]
   If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in
 …
 \xmlline|<context src="./nemo_def.xml" />|
+\noindent In NEMO, by default, the field and domain definition is done in 2 separate files:
+\path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that
+\noindent In NEMO, by default, the field definition is done in 3 separate files (
+\path{cfgs/SHARED/field_def_nemo-oce.xml},
+\path{cfgs/SHARED/field_def_nemo-pisces.xml} and
+\path{cfgs/SHARED/field_def_nemo-ice.xml} ) and the  domain definition is done in another file ( \path{cfgs/SHARED/domain_def_nemo.xml} )
+that
 are included in the main iodef.xml file through the following commands:
 \begin{xmllines}
+<field_definition src="./field_def.xml" />
+<domain_definition src="./domain_def.xml" />
+<context id="nemo" src="./context_nemo.xml"/>
 \end{xmllines}
 …
 Secondly, the group can be used to replace a list of elements.
+Several examples of groups of fields are proposed at the end of the file \path{CONFIG/SHARED/field_def.xml}.
+Several examples of groups of fields are proposed at the end of the XML field files (
+\path{cfgs/SHARED/field_def_nemo-oce.xml},
+\path{cfgs/SHARED/field_def_nemo-pisces.xml} and
+\path{cfgs/SHARED/field_def_nemo-ice.xml} ) .
 For example, a short list of the usual variables related to the U grid:
 …
 <field_group id="groupU" >
    <field field_ref="uoce"  />
    <field field_ref="suoce" />
+   <field field_ref="ssu" />
    <field field_ref="utau"  />
 </field_group>
 …
 the tag family domain.
 It must therefore be done in the domain part of the XML file.
 For example, in \path{CONFIG/SHARED/domain_def.xml}, we provide the following example of a definition of
+For example, in \path{cfgs/SHARED/domain_def.xml}, we provide the following example of a definition of
 a 5 by 5 box with the bottom left corner at point (10,10).
 \begin{xmllines}
 <domain_group id="grid_T">
    <domain id="myzoom" zoom_ibegin="10" zoom_jbegin="10" zoom_ni="5" zoom_nj="5" />
+<domain id="myzoomT" domain_ref="grid_T">
+   <zoom_domain ibegin="10" jbegin="10" ni="5" nj="5" />
 \end{xmllines}
 …
 \begin{xmllines}
 <file id="myfile_vzoom" output_freq="1d" >
    <field field_ref="toce" domain_ref="myzoom"/>
+   <field field_ref="toce" domain_ref="myzoomT"/>
 </file>
 \end{xmllines}
 …
 \subsubsection{Define vertical zooms}
 Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis.
+Vertical zooms are defined through the attributs zoom\_begin and zoom\_n of the tag family axis.
 It must therefore be done in the axis part of the XML file.
+For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example:
+\begin{xmllines}
+<axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" >
+   <axis id="deptht" />
+   <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" />
+For example, in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef_demo.xml}, we provide the following example:
+\begin{xmllines}
+<axis_definition>
+   <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" />
+   <axis id="deptht_zoom" azix_ref="deptht" >
+      <zoom_axis zoom_begin="1" zoom_n="10" />
+   </axis>
 \end{xmllines}
 …
 \end{xmllines}
 Note that, then the code is crashing, writting real4 variables forces a numerical convection from
+Note that, then the code is crashing, writting real4 variables forces a numerical conversion from
 real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files.
 Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem.
 …
     \hline
   \end{tabularx}
   \caption{Field tags ("\tt{field\_*}")}
+  \caption{Field tags ("\ttfamily{field\_*}")}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{File tags ("\tt{file\_*}")}
+  \caption{File tags ("\ttfamily{file\_*}")}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{Axis tags ("\tt{axis\_*}")}
+  \caption{Axis tags ("\ttfamily{axis\_*}")}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{Domain tags ("\tt{domain\_*)}"}
+  \caption{Domain tags ("\ttfamily{domain\_*)}"}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{Grid tags ("\tt{grid\_*}")}
+  \caption{Grid tags ("\ttfamily{grid\_*}")}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{Reference attributes ("\tt{*\_ref}")}
+  \caption{Reference attributes ("\ttfamily{*\_ref}")}
 \end{table}
 …
     \hline
   \end{tabularx}
   \caption{Domain attributes ("\tt{zoom\_*}")}
+  \caption{Domain attributes ("\ttfamily{zoom\_*}")}
 \end{table}
 …
 \subsection{CF metadata standard compliance}
 Output from the XIOS-1.0 IO server is compliant with
+Output from the XIOS IO server is compliant with
 \href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of
 the CF metadata standard.
 …
 %       NetCDF4 support
 % ================================================================
+\section{NetCDF4 support (\protect\key{netcdf4})}
+\section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})]
+{NetCDF4 support (\protect\key{netcdf4})}
 \label{sec:DIA_nc4}
 …
 Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead.
 For a fuller discussion on chunking and other performance issues the reader is referred to
 the NetCDF4 documentation found \href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}.
+the NetCDF4 documentation found \href{https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_perf_chunking.html}{here}.
 The new features are only available when the code has been linked with a NetCDF4 library
 …
 \end{forlines}
 \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in
 the mono-processor case (\ie global domain of {\small\tt 182x149x31}).
+\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\ttfamily 46x38x1} respectively in
+the mono-processor case (\ie global domain of {\small\ttfamily 182x149x31}).
 An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in
 table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with
 …
 %       Tracer/Dynamics Trends
 % -------------------------------------------------------------------------------------------------------------
+\section{Tracer/Dynamics trends  (\protect\ngn{namtrd})}
+\section[Tracer/Dynamics trends (\texttt{namtrd})]
+{Tracer/Dynamics trends (\protect\ngn{namtrd})}
 \label{sec:DIA_trd}
 …
 (\ie at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines).
 This capability is controlled by options offered in \ngn{namtrd} namelist.
 Note that the output are done with xIOS, and therefore the \key{IOM} is required.
+Note that the output are done with XIOS, and therefore the \key{iomput} is required.
 What is done depends on the \ngn{namtrd} logical set to \forcode{.true.}:
 …
 Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via
 the \key{trdtrc} and \key{trdmld\_trc} CPP keys.
+the \key{trdtrc} and \key{trdmxl\_trc} CPP keys.
 \textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested.
 In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working,
 and none of the options have been tested with variable volume (\ie \key{vvl} defined).
+and none of the options have been tested with variable volume (\ie \np{ln\_linssh}\forcode{ = .true.}).
 % -------------------------------------------------------------------------------------------------------------
 %       On-line Floats trajectories
 % -------------------------------------------------------------------------------------------------------------
+\section{FLO: On-Line Floats trajectories (\protect\key{floats})}
+\section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})]
+{FLO: On-Line Floats trajectories (\protect\key{floats})}
 \label{sec:FLO}
 %--------------------------------------------namflo-------------------------------------------------------
 …
 The on-line computation of floats advected either by the three dimensional velocity field or constraint to
 remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project.
 Options are defined by \ngn{namflo} namelis variables.
 The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option),
+Options are defined by \ngn{namflo} namelist variables.
+The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option),
 or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}).
 Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which
+Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which
 are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry.
 …
 In case of Ariane convention, input filename is \np{init\_float\_ariane}.
 Its format is: \\
 {\scriptsize \texttt{I J K nisobfl itrash itrash}}
+{\scriptsize \texttt{I J K nisobfl itrash}}
 \noindent with:
 …
 In that case, output filename is trajec\_float.
+Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}).
+There are 2 possibilities:
+- if (\key{iomput}) is used, outputs are selected in  iodef.xml.
+Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}) with
+\key{iomput} and outputs selected in iodef.xml.
 Here it is an example of specification to put in files description section:
 …
 \end{xmllines}
- -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library.
- See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of
- this marvellous diagnostic tool.
 % -------------------------------------------------------------------------------------------------------------
 %       Harmonic analysis of tidal constituents
 % -------------------------------------------------------------------------------------------------------------
+\section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) }
+\section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})]
+{Harmonic analysis of tidal constituents (\protect\key{diaharm})}
 \label{sec:DIA_diag_harm}
 %------------------------------------------namdia_harm----------------------------------------------------
+%------------------------------------------nam_diaharm----------------------------------------------------
+%
 \nlst{nam_diaharm}
 …
 This on-line Harmonic analysis is actived with \key{diaharm}.
 Some parameters are available in namelist \ngn{namdia\_harm}:
+Some parameters are available in namelist \ngn{nam\_diaharm}:
  - \np{nit000\_han} is the first time step used for harmonic analysis
 …
 %       Sections transports
 % -------------------------------------------------------------------------------------------------------------
+\section{Transports across sections (\protect\key{diadct}) }
+\section[Transports across sections (\texttt{\textbf{key\_diadct}})]
+{Transports across sections (\protect\key{diadct})}
 \label{sec:DIA_diag_dct}
 …
 Each section is defined by the coordinates of its 2 extremities.
 The pathways between them are contructed using tools which can be found in \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT}
 and are written in a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is later read in by
+The pathways between them are contructed using tools which can be found in \texttt{tools/SECTIONS\_DIADCT}
+and are written in a binary file \texttt{section\_ijglobal.diadct} which is later read in by
 NEMO to compute on-line transports.
 …
 \subsubsection{Creating a binary file containing the pathway of each section}
 In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run},
+In \texttt{tools/SECTIONS\_DIADCT/run},
 the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed
 (this list of sections is based on MERSEA project metrics).
 …
  The script \texttt{job.ksh} computes the pathway for each section and creates a binary file
  \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\
+ \texttt{section\_ijglobal.diadct} which is read by NEMO. \\
  It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with
 …
 The steric effect is therefore not explicitely represented.
 This approximation does not represent a serious error with respect to the flow field calculated by the model
 \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level,
+\citep{greatbatch_JGR94}, but extra attention is required when investigating sea level,
 as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales.
 This is especially true for investigation into sea level rise due to global warming.
 Fortunately, the steric contribution to the sea level consists of a spatially uniform component that
 can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}.
+can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}.
 In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed,
 we compare, in the following, the non-Boussinesq and Boussinesq cases.
 …
 the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid.
 Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by
+Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by
 considering the mass budget of the ocean.
 The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux
 must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean
 \citep{Greatbatch_JGR94}.
+\citep{greatbatch_JGR94}.
 In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$,
 the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level,
 …
 This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since,
 with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than
 $\%$ from this value (\cite{Gill1982}, page 47).
+$\%$ from this value (\cite{gill_bk82}, page 47).
 Second, we have assumed here that the total ocean surface, $\mathcal{A}$,
 …
 Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered.
 In the non linear free surface case, \ie \key{vvl} defined, it is given by
+In the non linear free surface case, \ie \np{ln\_linssh}\forcode{ = .true.}, it is given by
 \[
 …
 so that there are no associated ocean currents.
 Hence, the dynamically relevant sea level is the effective sea level,
 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}.
+\ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}.
 However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between
 ice and ocean.
 …
 where $S_o$ and $p_o$ are the initial salinity and pressure, respectively.
+Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to
+be called.
+Both steric and thermosteric sea level are computed in \mdl{diaar5}.
 % -------------------------------------------------------------------------------------------------------------
 %       Other Diagnostics
 % -------------------------------------------------------------------------------------------------------------
+\section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})}
+\section[Other diagnostics]
+{Other diagnostics}
 \label{sec:DIA_diag_others}
 …
 The available ready-to-add diagnostics modules can be found in directory DIA.
+\subsection{Depth of various quantities (\protect\mdl{diahth})}
+\subsection[Depth of various quantities (\textit{diahth.F90})]
+{Depth of various quantities (\protect\mdl{diahth})}
 Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key:
 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
+- the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth})
 - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
 …
 - the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
-% -----------------------------------------------------------
-%     Poleward heat and salt transports
-% -----------------------------------------------------------
-\subsection{Poleward heat and salt transports (\protect\mdl{diaptr})}
-%------------------------------------------namptr-----------------------------------------
-\nlst{namptr}
-%-----------------------------------------------------------------------------------------
-The poleward heat and salt transports, their advective and diffusive component,
-and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true
-(see the \textit{\ngn{namptr} } namelist below).
-When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian,
-Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean.
-The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,
-the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_mask_subasins}
+    \includegraphics[width=\textwidth]{Fig_mask_subasins}
     \caption{
       \protect\label{fig:mask_subasins}
 …
 %       CMIP specific diagnostics
 % -----------------------------------------------------------
+\subsection{CMIP specific diagnostics (\protect\mdl{diaar5})}
+A series of diagnostics has been added in the \mdl{diaar5}.
+They corresponds to outputs that are required for AR5 simulations (CMIP5)
+\subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})]
+{CMIP specific diagnostics (\protect\mdl{diaar5})}
+A series of diagnostics has been added in the \mdl{diaar5} and \mdl{diaptr}.
+In \mdl{diaar5} they correspond to outputs that are required for AR5 simulations (CMIP5)
 (see also \autoref{sec:DIA_steric} for one of them).
+Activating those outputs requires to define the \key{diaar5} CPP key.
+The module \mdl{diaar5} is active when one of the following outputs is required : global total volume (voltot), global mean ssh (sshtot), global total mass (masstot), global mean temperature (temptot), global mean ssh steric (sshsteric), global mean ssh thermosteric (sshthster), global mean salinity (saltot), sea water pressure at sea floor (botpres), dynamic sea surface height (sshdyn).
+In \mdl{diaptr} when \np{ln\_diaptr}\forcode{ = .true.}
+(see the \textit{\ngn{namptr} } namelist below) can be computed on-line the poleward heat and salt transports, their advective and diffusive component, and the meriodional stream function .
+When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian,
+Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean.
+The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays,
+the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}).
+%------------------------------------------namptr-----------------------------------------
+\nlst{namptr}
+%-----------------------------------------------------------------------------------------
 % -----------------------------------------------------------

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIU.tex

-                      r10442
+                      r11422
 (\ie from the temperature of the top few model levels) or from some other source.
 It must be noted that both the cool skin and warm layer models produce estimates of the change in temperature
 ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) and
+($\Delta T_{\mathrm{cs}}$ and $\Delta T_{\mathrm{wl}}$) and
 both must be added to a foundation SST to obtain the true skin temperature.
 …
 %===============================================================
 The warm layer is calculated using the model of \citet{Takaya_al_JGR10} (TAKAYA10 model hereafter).
+The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter).
 This is a simple flux based model that is defined by the equations
 \begin{align}
 \frac{\partial{\Delta T_{\rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p
+\frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p
 \nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,}
 \label{eq:ecmwf1} \\
 L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2}
 \end{align}
 where $\Delta T_{\rm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal.
+where $\Delta T_{\mathrm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal.
 In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water,
 $\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water,
 $\rho_w$ is the water density, and $L$ is the Monin-Obukhov length.
 The tunable variable $\nu$ is a shape parameter that defines the expected subskin temperature profile via
 $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\rm{wl}}$,
+$T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\mathrm{wl}}$,
 where $T$ is the absolute temperature and $z\le D_T$ is the depth below the top of the warm layer.
 The influence of wind on TAKAYA10 comes through the magnitude of the friction velocity of the water $u^*_{w}$,
 …
 the diurnal layer, \ie
 \[
   Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,}
+  Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,}
   % \label{eq:e_flux_eqn}
 \]
 where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux,
 and $Q_{\rm{sol}}$ is the solar flux absorbed within the diurnal warm layer.
 For $Q_{\rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used.
+where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux,
+and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer.
+For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used.
 In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$,
 where $L_a=0.3$\footnote{
 …
 %===============================================================
 The cool skin is modelled using the framework of \citet{Saunders_JAS82} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$.
 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\rm{cs}}$ becomes
+The cool skin is modelled using the framework of \citet{saunders_JAS67} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$.
+As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\mathrm{cs}}$ becomes
 \[
   % \label{eq:sunders_eqn}
   \Delta T_{\rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,}
+  \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,}
 \]
 where $Q_{\rm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and
+where $Q_{\mathrm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and
 $k_t$ is the thermal conductivity of sea water.
 $\delta$ is the thickness of the skin layer and is given by
 …
 \end{equation}
 where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which
 \citet{Saunders_JAS82} suggested varied between 5 and 10.
+\citet{saunders_JAS67} suggested varied between 5 and 10.
 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{Artale_al_JGR02},
 which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other parametrisations at
+The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02},
+which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at
 both low and high wind speeds.
 Specifically,

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DOM.tex

-                      r10502
+                      r11422
 %     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled
+\vfill
+\begin{figure}[b]
+\subsubsection*{Changes record}
+\begin{tabular}{m{0.08\linewidth}||m{0.32\linewidth}|m{0.6\linewidth}}
+    Release   & Author(s)     & Modifications \\
+\hline
+    {\em 4.0} & {\em Simon M{\"u}ller \& Andrew Coward} & {\em Compatibility changes for v4.0. Major simplication has moved many of the options to external domain configuration tools. For now this information has been retained in \autoref{apdx:DOMAINcfg} }  \\
+    {\em 3.x} & {\em Sebastien Masson, Gurvan Madec \& Rashid Benshila } & {\em }  \\
+\end{tabular}
+\end{figure}
 \newpage
 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP},
 we need to choose a discretization on a grid, and numerical algorithms.
+we need to choose a grid for spatial discretization and related numerical algorithms.
 In the present chapter, we provide a general description of the staggered grid used in \NEMO,
 and other information relevant to the main directory routines as well as the DOM (DOMain) directory.
+and other relevant information about the DOM (DOMain) source-code modules .
 % ================================================================
 …
 \begin{figure}[!tb]
   \begin{center}
     \includegraphics[]{Fig_cell}
+    \includegraphics[width=\textwidth]{Fig_cell}
     \caption{
       \protect\label{fig:cell}
 …
 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional,
 centred second-order finite difference approximation.
 Special attention has been given to the homogeneity of the solution in the three space directions.
+Special attention has been given to the homogeneity of the solution in the three spatial directions.
 The arrangement of variables is the same in all directions.
 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in
 the centre of each face of the cells (\autoref{fig:cell}).
 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification
 \citep{Mesinger_Arakawa_Bk76}.
+\citep{mesinger.arakawa_bk76}.
 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and
 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.
 …
 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}.
 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
 $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity.
+$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity.
 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation
 while the scale factors are chosen equal to their local analytical value.
 …
 the continuous properties (see \autoref{apdx:C}).
 A similar, related remark can be made about the domain size:
 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors
+when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors
 (see \autoref{eq:DOM_bar} in the next section).
 …
     \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
       \hline
       T  & $i      $ & $j      $ & $k      $ \\
+      t  & $i      $ & $j      $ & $k      $ \\
       \hline
       u  & $i + 1/2$ & $j      $ & $k      $ \\
 …
       \protect\label{tab:cell}
       Location of grid-points as a function of integer or integer and a half value of the column, line or level.
       This indexing is only used for the writing of the semi -discrete equation.
       In the code, the indexing uses integer values only and has a reverse direction in the vertical
+      This indexing is only used for the writing of the semi -discrete equations.
+      In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface.
       (see \autoref{subsec:DOM_Num_Index})
+    }
   \end{center}
 \end{table}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+Note that the definition of the scale factors
+(\ie as the analytical first derivative of the transformation that
+results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
+is specific to the \NEMO model \citep{marti.madec.ea_JGR92}.
+As an example, a scale factor in the $i$ direction is defined locally at a $t$-point,
+whereas many other models on a C grid choose to define such a scale factor as
+the distance between the $u$-points on each side of the $t$-point.
+Relying on an analytical transformation has two advantages:
+firstly, there is no ambiguity in the scale factors appearing in the discrete equations,
+since they are first introduced in the continuous equations;
+secondly, analytical transformations encourage good practice by the definition of smoothly varying grids
+(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}.
+An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[width=\textwidth]{Fig_zgr_e3}
+    \caption{
+      \protect\label{fig:zgr_e3}
+      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
+      and (b) analytically derived grid-point position and scale factors.
+      For both grids here, the same $w$-point depth has been chosen but
+      in (a) the $t$-points are set half way between $w$-points while
+      in (b) they are defined from an analytical function:
+      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
+      Note the resulting difference between the value of the grid-size $\Delta_k$ and
+      those of the scale factor $e_k$.
+    }
+  \end{center}
+\end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at
 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while
 its Laplacian is defined at $t$-point.
+its Laplacian is defined at the $t$-point.
 These operators have the following discrete forms in the curvilinear $s$-coordinates system:
 \[
 …
 \end{equation}
 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which
 is a masked field (i.e. equal to zero inside solid area):
+The vertical average over the whole water column is denoted by an overbar and is for
+a masked field $q$ (\ie a quantity that is equal to zero inside solid areas):
 \begin{equation}
   \label{eq:DOM_bar}
 …
 \end{equation}
 where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points,
 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over
+$k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over
 all grid points of the same type in the direction indicated by the subscript (here $k$).
 …
 vector points $(u,v,w)$.
 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area.
 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$)
+Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas.
+It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$)
 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$,
 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie
 …
 \begin{figure}[!tb]
   \begin{center}
     \includegraphics[]{Fig_index_hor}
+    \includegraphics[width=\textwidth]{Fig_index_hor}
     \caption{
       \protect\label{fig:index_hor}
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 The array representation used in the \fortran code requires an integer indexing while
 the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of
 integer values for $t$-points and both integer and integer and a half values for all the other points.
 Therefore a specific integer indexing must be defined for points other than $t$-points
+The array representation used in the \fortran code requires an integer indexing.
+However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of
+integer values for $t$-points only while all the other points involve integer and a half values.
+Therefore, a specific integer indexing has been defined for points other than $t$-points
 (\ie velocity and vorticity grid-points).
 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$.
+Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$.
 % -----------------------------------
 …
 \label{subsec:DOM_Num_Index_vertical}
 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in
 the \fortran code compared to the indexing used in the semi -discrete equations and
+In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in
+the \fortran code is the reverse of that used in the semi -discrete equations and
 given in \autoref{subsec:DOM_cell}.
 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below
+The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below
 (\autoref{fig:index_vert}).
+The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while
+the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}).
+Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index,
+in opposition to what is done in the horizontal plane where
+it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that
+have the same $i$ or $j$ index
+The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while
+the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}).
+Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie $t$-points and their
+nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where
+the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index
 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}).
 Since the scale factors are chosen to be strictly positive,
 a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of
 the discrete equations given in this documentation.
+a \textit{minus sign} is included in the \fortran implementations of \textit{all the vertical derivatives} of
+the discrete equations given in this manual in order to accommodate the opposing vertical index directions in implementation and documentation.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!pt]
   \begin{center}
     \includegraphics[]{Fig_index_vert}
+    \includegraphics[width=\textwidth]{Fig_index_vert}
     \caption{
       \protect\label{fig:index_vert}
       Vertical integer indexing used in the \fortran code.
       Note that the $k$-axis is orientated downward.
       The dashed area indicates the cell in which variables contained in arrays have the same $k$-index.
+      Note that the $k$-axis is oriented downward.
+      The dashed area indicates the cell in which variables contained in arrays have a common $k$-index.
+    }
   \end{center}
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% -------------------------------------------------------------------------------------------------------------
+%        Domain configuration
+% -------------------------------------------------------------------------------------------------------------
+\section{Spatial domain configuration}
+\label{subsec:DOM_config}
+\nlst{namcfg}
+Two typical methods are available to specify the spatial domain
+configuration; they can be selected using parameter \np{ln\_read\_cfg}
+parameter in namelist \ngn{namcfg}.
+If \np{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters
+and fields are read from a netCDF input file, whose name (without its .nc
+suffix) can be specified as the value of the \np{cn\_domcfg} parameter in
+namelist \ngn{namcfg}.
+If \np{ln\_read\_cfg} is set to \forcode{.false.}, the domain-specific
+parameters and fields can be provided (\eg analytically computed) by subroutines
+\mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in
+the \path{MY_SRC} directory of the configuration, and default versions that
+configure the spatial domain for the GYRE reference configuration are present in
+the \path{src/OCE/USR} directory.
+In version 4.0 there are no longer any options for reading complex bathmetries and
+performing a vertical discretization at run-time. Whilst it is occasionally convenient
+to have a common bathymetry file and, for example, to run similar models with and
+without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to
+overly complex code. Worse still is the difficulty of ensuring the model configurations
+intended to be identical are indeed so when the model domain itself can be altered by runtime
+selections. The code previously used to perform vertical discretization has be incorporated
+into an external tool (\path{tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}.
+The next subsections summarise the parameter and fields related to the
+configuration of the whole model domain. These represent the minimum information
+that must be provided either via the \np{cn\_domcfg} file or set by code
+inserted into user-supplied versions of the \mdl{usrdef\_*} subroutines. The
+requirements are presented in three sections: the domain size
+(\autoref{subsec:DOM_size}), the horizontal mesh
+(\autoref{subsec:DOM_hgr}), and the vertical grid
+(\autoref{subsec:DOM_zgr}).
 % -----------------------------------
 %        Domain Size
 % -----------------------------------
 \subsubsection{Domain size}
+\subsection{Domain size}
 \label{subsec:DOM_size}
+The total size of the computational domain is set by the parameters \np{jpiglo},
+\np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively.
+Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when
+the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined,
+see \autoref{sec:LBC_mpp}).
+% ================================================================
+% Domain: List of fields needed
+% ================================================================
+\section{Needed fields}
+\label{sec:DOM_fields}
+The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that
+gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
+The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}.
+The associated scale factors are defined using the analytical first derivative of the transformation
+\autoref{eq:scale_factors}.
+Necessary fields for configuration definition are:
+\begin{itemize}
+\item
+  Geographic position:
+  longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and
+  latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif}
+  (all respectively at T, U, V and F point)
+\item
+  Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t}
+  (at T and F point)
+\item
+  Scale factors:
+  \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction),
+  \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and
+  \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\
+  \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits),
+  \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed.
+\end{itemize}
+These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in
+\ngn{namcfg}.
+\nlst{namcfg}
+Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration.
+For Reference Configurations of NEMO input domain files are supplied by NEMO System Team.
+For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}.
+They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory,
+they provide the horizontal and vertical mesh.
+% -------------------------------------------------------------------------------------------------------------
+%        Needed fields
+% -------------------------------------------------------------------------------------------------------------
+%\subsection{List of needed fields to build DOMAIN}
+%\label{subsec:DOM_fields_list}
+The total size of the computational domain is set by the parameters
+\np{jpiglo}, \np{jpjglo} and \np{jpkglo} for the $i$, $j$ and $k$
+directions, respectively. Note, that the variables \forcode{jpi} and \forcode{jpj}
+refer to the size of each processor subdomain when the code is run in
+parallel using domain decomposition (\key{mpp\_mpi} defined, see
+\autoref{sec:LBC_mpp}).
+The name of the configuration is set through parameter \np{cn\_cfg},
+and the nominal resolution through parameter \np{nn\_cfg} (unless in
+the input file both of variables \forcode{ORCA} and \forcode{ORCA_index}
+are present, in which case \np{cn\_cfg} and \np{nn\_cfg} are set from these
+values accordingly).
+The global lateral boundary condition type is selected from 8 options
+using parameter \np{jperio}. See \autoref{sec:LBC_jperio} for
+details on the available options and the corresponding values for
+\np{jperio}.
 % ================================================================
 % Domain: Horizontal Grid (mesh)
 % ================================================================
+\section{Horizontal grid mesh (\protect\mdl{domhgr})}
+\label{sec:DOM_hgr}
+% -------------------------------------------------------------------------------------------------------------
+%        Coordinates and scale factors
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Coordinates and scale factors}
+\label{subsec:DOM_hgr_coord_e}
+The ocean mesh (\ie the position of all the scalar and vector points) is defined by
+the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
+The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}.
+The associated scale factors are defined using the analytical first derivative of the transformation
+\autoref{eq:scale_factors}.
+These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr},
+which provide the horizontal and vertical meshes, respectively.
+This section deals with the horizontal mesh parameters.
+In a horizontal plane, the location of all the model grid points is defined from
+the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$.
+The horizontal scale factors are calculated using \autoref{eq:scale_factors}.
+For example, when the longitude and latitude are function of a single value
+($i$ and $j$, respectively) (geographical configuration of the mesh),
+the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$,
+and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module.
+The model computes the grid-point positions and scale factors in the horizontal plane as follows:
+\begin{align*}
+   \lambda_t &\equiv \text{glamt} =      \lambda (i      )
+  &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\
+   \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2)
+  &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\
+   \lambda_v &\equiv \text{glamv} =      \lambda (i      )
+  &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\
+   \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2)
+  &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\
+   e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) |
+  &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
+   e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) |
+  &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
+   e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) |
+  &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\
+   e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) |
+  &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         |
+\end{align*}
+where the last letter of each computational name indicates the grid point considered and
+$r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants).
+Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points,
+thus no specific arrays are defined at $w$-points.
+Note that the definition of the scale factors
+(\ie as the analytical first derivative of the transformation that
+gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
+is specific to the \NEMO model \citep{Marti_al_JGR92}.
+As an example, $e_{1t}$ is defined locally at a $t$-point,
+whereas many other models on a C grid choose to define such a scale factor as
+the distance between the $U$-points on each side of the $t$-point.
+Relying on an analytical transformation has two advantages:
+firstly, there is no ambiguity in the scale factors appearing in the discrete equations,
+since they are first introduced in the continuous equations;
+secondly, analytical transformations encourage good practice by the definition of smoothly varying grids
+(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}.
+An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]
+  \begin{center}
+    \includegraphics[]{Fig_zgr_e3}
+    \caption{
+      \protect\label{fig:zgr_e3}
+      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
+      and (b) analytically derived grid-point position and scale factors.
+      For both grids here, the same $w$-point depth has been chosen but
+      in (a) the $t$-points are set half way between $w$-points while
+      in (b) they are defined from an analytical function:
+      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
+      Note the resulting difference between the value of the grid-size $\Delta_k$ and
+      those of the scale factor $e_k$.
+    }
+  \end{center}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% -------------------------------------------------------------------------------------------------------------
+%        Choice of horizontal grid
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Choice of horizontal grid}
+\label{subsec:DOM_hgr_msh_choice}
+% -------------------------------------------------------------------------------------------------------------
+%        Grid files
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Output grid files}
+\label{subsec:DOM_hgr_files}
+All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks)
+can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}).
+This can be particularly useful for plots and off-line diagnostics.
+In some cases, the user may choose to make a local modification of a scale factor in the code.
+This is the case in global configurations when restricting the width of a specific strait
+(usually a one-grid-point strait that happens to be too wide due to insufficient model resolution).
+An example is Gibraltar Strait in the ORCA2 configuration.
+When such modifications are done,
+the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid.
+\subsection{Horizontal grid mesh (\protect\mdl{domhgr})}
+\label{subsec:DOM_hgr}
+% ================================================================
+% Domain: List of hgr-related fields needed
+% ================================================================
+\subsubsection{Required fields}
+\label{sec:DOM_hgr_fields}
+The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include:
+\begin{Verbatim}[fontsize=\tiny]
+int    jpiglo, jpjglo, jpkglo            /* global domain sizes                                          */
+int    jperio                            /* lateral global domain b.c.                                   */
+double glamt, glamu, glamv, glamf        /* geographic longitude (t,u,v and f points respectively)       */
+double gphit, gphiu, gphiv, gphif        /* geographic latitude                                          */
+double e1t, e1u, e1v, e1f                /* horizontal scale factors                                     */
+double e2t, e2u, e2v, e2f                /* horizontal scale factors                                     */
+\end{Verbatim}
+The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, evaluated at the values as specified in Table \autoref{tab:cell} for the respective grid-point position. The calculation of the values of the horizontal scale factor arrays in general additionally involves partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, evaluated for the same arguments as $\lambda$ and $\varphi$.
+\subsubsection{Optional fields}
+\begin{Verbatim}[fontsize=\tiny]
+                                         /* Optional:                                                    */
+int    ORCA, ORCA_index                  /* configuration name, configuration resolution                 */
+double e1e2u, e1e2v                      /* U and V surfaces (if grid size reduction in some straits)    */
+double ff_f, ff_t                        /* Coriolis parameter (if not on the sphere)                    */
+\end{Verbatim}
+NEMO can support the local reduction of key strait widths by altering individual values of
+e2u or e1v at the appropriate locations. This is particularly useful for locations such as
+Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for
+illustrated examples). The key is to reduce the faces of $T$-cell (\ie change the value of
+the horizontal scale factors at $u$- or $v$-point) but not the volume of the cells. Doing
+otherwise can lead to numerical instability issues.  In normal operation the surface areas
+are computed from $\texttt{e1u} * \texttt{e2u}$ and $\texttt{e1v} * \texttt{e2v}$ but in
+cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$
+grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed
+in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn\_domcfg} file they are
+read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set
+their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area
+computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation.
+\smallskip
+Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t}
+which can be used to provide the Coriolis parameter at F- and T-points respectively if the
+mesh is not on a sphere. If present these fields will be read and used and the normal
+calculation ($2*\Omega*\sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set
+their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation
+flag: \texttt{iff} to a non-zero value to suppress their re-computation.
+Note that longitudes, latitudes, and scale factors at $w$ points are exactly
+equal to those of $t$ points, thus no specific arrays are defined at $w$ points.
 % ================================================================
 % Domain: Vertical Grid (domzgr)
 % ================================================================
+\section{Vertical grid (\protect\mdl{domzgr})}
+\label{sec:DOM_zgr}
+%-----------------------------------------nam_zgr & namdom-------------------------------------------
+%
+%\nlst{namzgr}
+\nlst{namdom}
+\subsection[Vertical grid (\textit{domzgr.F90})]
+{Vertical grid (\protect\mdl{domzgr})}
+\label{subsec:DOM_zgr}
+%-----------------------------------------namdom-------------------------------------------
+\nlst{namdom}
 %-------------------------------------------------------------------------------------------------------------
-Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists.
 In the vertical, the model mesh is determined by four things:
+(1) the bathymetry given in meters;
+(2) the number of levels of the model (\jp{jpk});
+(3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and
+(4) the masking system, \ie the number of wet model levels at each
+$(i,j)$ column of points.
+\begin{enumerate}
+  \item the bathymetry given in meters;
+  \item the number of levels of the model (\jp{jpk});
+  \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and
+  \item the masking system, \ie the number of wet model levels at each
+$(i,j)$ location of the horizontal grid.
+\end{enumerate}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb]
   \begin{center}
     \includegraphics[]{Fig_z_zps_s_sps}
+    \includegraphics[width=\textwidth]{Fig_z_zps_s_sps}
     \caption{
       \protect\label{fig:z_zps_s_sps}
 …
       (d) hybrid $s-z$ coordinate,
       (e) hybrid $s-z$ coordinate with partial step, and
       (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}).
+      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}).
       Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).
+    }
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,
+must be done once of all at the beginning of an experiment.
+It is not intended as an option which can be enabled or disabled in the middle of an experiment.
+Three main choices are offered (\autoref{fig:z_zps_s_sps}):
+$z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}),
+$z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}),
+or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}).
+Hybridation of the three main coordinates are available:
+$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}).
+By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that
+the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}).
+When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}),
+the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface
+(in other words, the top of the ocean in not a rigid-lid).
+The last choice in terms of vertical coordinate concerns the presence (or not) in
+the model domain of ocean cavities beneath ice shelves.
+Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in.
+This option is currently only available in $z$- or $zps$-coordinate,
+and partial step are also applied at the ocean/ice shelf interface.
+Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for
+reading it from a file.
+The only input file is the bathymetry (in meters) (\ifile{bathy\_meter})
+\footnote{
+  N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file,
+  so that the computation of the number of wet ocean point in each water column is by-passed}.
+If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing
+the ice shelf draft (in meters) is needed.
+After reading the bathymetry, the algorithm for vertical grid definition differs between the different options:
+\begin{description}
+\item[\textit{zco}]
+  set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$.
+\item[\textit{zps}]
+  set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at
+  each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
+\item[\textit{sco}]
+  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
+  set the three-dimensional transformation.
+\item[\textit{s-z} and \textit{s-zps}]
+  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
+  set the three-dimensional transformation $z(i,j,k)$,
+  and possibly introduce masking of extra land points to better fit the original bathymetry file.
+\end{description}
+%%%
+\gmcomment{   add the description of the smoothing:  envelop topography...}
+%%%
+Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}),
+the arrays describing the grid point depths and vertical scale factors are three set of
+three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step.
+The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.
+They are updated at each model time step using a fixed reference coordinate system which
+computer names have a $\_0$ suffix.
+When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before},
+\textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart.
+% -------------------------------------------------------------------------------------------------------------
+%        Meter Bathymetry
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Meter bathymetry}
+\label{subsec:DOM_bathy}
+Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy}
+(found in \ngn{namdom} namelist):
+\begin{description}
+\item[\np{nn\_bathy}~\forcode{= 0}]:
+  a flat-bottom domain is defined.
+  The total depth $z_w (jpk)$ is given by the coordinate transformation.
+  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}.
+\item[\np{nn\_bathy}~\forcode{= -1}]:
+  a domain with a bump of topography one third of the domain width at the central latitude.
+  This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount.
+\item[\np{nn\_bathy}~\forcode{= 1}]:
+  read a bathymetry and ice shelf draft (if needed).
+  The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at
+  each grid point of the model grid.
+  The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto
+  the horizontal ocean mesh.
+  Defining the bathymetry also defines the coastline: where the bathymetry is zero,
+  no model levels are defined (all levels are masked).
+  The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at
+  each grid point of the model grid.
+  This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}.
+  Defining the ice shelf draft will also define the ice shelf edge and the grounding line position.
+\end{description}
+When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies
+(\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with
+the rest of the ocean.
+This is unnecessary when the ocean is forced by fixed atmospheric conditions,
+so these seas can be removed from the ocean domain.
+The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}),
+but the code has to be adapted to the user's configuration.
+% -------------------------------------------------------------------------------------------------------------
+%        z-coordinate  and reference coordinate transformation
+% -------------------------------------------------------------------------------------------------------------
+\subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate]
+            {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate}
+\label{subsec:DOM_zco}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!tb]
+  \begin{center}
+    \includegraphics[]{Fig_zgr}
+    \caption{
+      \protect\label{fig:zgr}
+      Default vertical mesh for ORCA2: 30 ocean levels (L30).
+      Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from
+      \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.
+    }
+  \end{center}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points,
+respectively.
+As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels.
+$gdepw_0(1)$ is the ocean surface.
+There are at most \jp{jpk}-1 $t$-points inside the ocean,
+the additional $t$-point at $jk = jpk$ is below the sea floor and is not used.
+The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose
+analytical derivative with respect to $k$ provides the vertical scale factors.
+The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$.
+This is done in routine \mdl{domzgr} through statement functions,
+using parameters provided in the \ngn{namcfg} namelist.
+It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}).
+In that case, the parameters \jp{jpk} (number of $w$-levels) and
+\np{pphmax} (total ocean depth in meters) fully define the grid.
+For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface.
+The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps):
+\begin{gather}
+  \label{eq:DOM_zgr_ana_1}
+    z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\
+    e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt|
+\end{gather}
+where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels.
+Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with
+a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}).
+If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same.
+However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to:
+\begin{equation}
+  \label{eq:DOM_zgr_ana_2}
+  \begin{split}
+    e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\
+    e_3^W(k) &= z_T (k    ) - z_T (k - 1)
+  \end{split}
+\end{equation}
+This formulation decrease the self-generated circulation into the ice shelf cavity
+(which can, in extreme case, leads to blow up).\\
+The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and
+a depth which varies from 0 at the sea surface to a minimum of $-5000~m$.
+This leads to the following conditions:
+\begin{equation}
+  \label{eq:DOM_zgr_coef}
+  \begin{array}{ll}
+    e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\
+    e_3 (jpk - 1/2) = 500. & z(jpk) = -5000.
+  \end{array}
+\end{equation}
+With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$,
+the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in
+\autoref{eq:DOM_zgr_ana_2} have been determined such that
+\autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method.
+For the first standard ORCA2 vertical grid this led to the following values:
+$h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$.
+The resulting depths and scale factors as a function of the model levels are shown in
+\autoref{fig:zgr} and given in \autoref{tab:orca_zgr}.
+Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.
+Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them.
+In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$.,
+in \ngn{namcfg} namelist, and specifies instead the four following parameters:
+The choice of a vertical coordinate is made when setting up the configuration;
+it is not intended to be an option which can be changed in the middle of an
+experiment. The one exception to this statement being the choice of linear or
+non-linear free surface. In v4.0 the linear free surface option is implemented
+as a special case of the non-linear free surface. This is computationally
+wasteful since it uses the structures for time-varying 3D metrics for fields
+that (in the linear free surface case) are fixed. However, the linear
+free-surface is rarely used and implementing it this way means a single configuration
+file can support both options.
+By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{ =
+.false.} in \ngn{namdom}): the coordinate follow the time-variation of the free
+surface so that the transformation is time dependent: $z(i,j,k,t)$
+(\eg \autoref{fig:z_zps_s_sps}f).  When a linear free surface is assumed
+(\np{ln\_linssh} set to \forcode{ = .true.} in \ngn{namdom}), the vertical
+coordinates are fixed in time, but the seawater can move up and down across the
+$z_0$ surface (in other words, the top of the ocean in not a rigid lid).
+Note that settings: \np{ln\_zco}, \np{ln\_zps}, \np{ln\_sco} and \np{ln\_isfcav} mentioned
+in the following sections appear to be namelist options but they are no longer truly
+namelist options for NEMO. Their value is written to and read from the domain configuration file
+and they should be treated as fixed parameters for a particular configuration. They are
+namelist options for the \forcode{DOMAINcfg} tool that can be used to build the
+configuration file and serve both to provide a record of the choices made whilst building the
+configuration and to trigger appropriate code blocks within NEMO.
+These values should not be altered in the \np{cn\_domcfg} file.
+\medskip
+The decision on these choices must be made when the \np{cn\_domcfg} file is constructed.
+Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c):
 \begin{itemize}
+\item
+  \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional).
+  The larger \np{ppacr}, the smaller the stretching.
+  Values from $3$ to $10$ are usual.
+\item
+  \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs
+  (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
+\item
+  \np{ppdzmin}: minimum thickness for the top layer (in meters).
+\item
+  \np{pphmax}: total depth of the ocean (meters).
+\item $z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}),
+\item $z$-coordinate with partial step ($zps$) bathymetry (\np{ln\_zps}\forcode{ = .true.}),
+\item Generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}).
 \end{itemize}
+As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are:
+\jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{table}
+  \begin{center}
+    \begin{tabular}{c||r|r|r|r}
+      \hline
+      \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\
+      \hline
+              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\
+      \hline
+              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\
+      \hline
+              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\
+      \hline
+              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\
+      \hline
+              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\
+      \hline
+              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\
+      \hline
+              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\
+      \hline
+             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\
+      \hline
+             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\
+      \hline
+             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\
+      \hline
+             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\
+      \hline
+             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\
+      \hline
+             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\
+      \hline
+             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\
+      \hline
+             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\
+      \hline
+             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\
+      \hline
+             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\
+      \hline
+             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\
+      \hline
+             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\
+      \hline
+             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\
+      \hline
+             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\
+      \hline
+             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\
+      \hline
+             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\
+      \hline
+             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\
+      \hline
+             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\
+      \hline
+             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\
+      \hline
+             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\
+      \hline
+             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\
+      \hline
+             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\
+      \hline
+    \end{tabular}
+  \end{center}
+  \caption{
+    \protect\label{tab:orca_zgr}
+    Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from
+    \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef}
+  }
+\end{table}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% -------------------------------------------------------------------------------------------------------------
+%        z-coordinate with partial step
+% -------------------------------------------------------------------------------------------------------------
+\subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})}
+\label{subsec:DOM_zps}
+%--------------------------------------------namdom-------------------------------------------------------
+\nlst{namdom}
+%--------------------------------------------------------------------------------------------------------------
+In $z$-coordinate partial step,
+the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in
+the previous section, \textit{except} in the bottom layer.
+The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to
+allow a better representation of the bathymetry, especially in the case of small slopes
+(where the bathymetry varies by less than one level thickness from one grid point to the next).
+The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry.
+With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$.
+The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$:
+the maximum thickness allowed is $2*e_{3t}(jpk - 1)$.
+This has to be kept in mind when specifying values in \ngn{namdom} namelist,
+as the maximum depth \np{pphmax} in partial steps:
+for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid,
+the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$).
+Two variables in the namdom namelist are used to define the partial step vertical grid.
+The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is
+the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat}
+(a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$).
+\gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  }
+% -------------------------------------------------------------------------------------------------------------
+%        s-coordinate
+% -------------------------------------------------------------------------------------------------------------
+\subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})}
+\label{subsec:DOM_sco}
+%------------------------------------------nam_zgr_sco---------------------------------------------------
+%
+%\nlst{namzgr_sco}
+%--------------------------------------------------------------------------------------------------------------
+Options are defined in \ngn{namzgr\_sco}.
+In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from
+the product of a depth field and either a stretching function or its derivative, respectively:
+\begin{align*}
+  % \label{eq:DOM_sco_ana}
+  z(k)   &= h(i,j) \; z_0 (k) \\
+  e_3(k) &= h(i,j) \; z_0'(k)
+\end{align*}
+where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and
+$z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom.
+The depth field $h$ is not necessary the ocean depth,
+since a mixed step-like and bottom-following representation of the topography can be used
+(\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}).
+The namelist parameter \np{rn\_rmax} determines the slope at which
+the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.
+The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as
+the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
+Options for stretching the coordinate are provided as examples,
+but care must be taken to ensure that the vertical stretch used is appropriate for the application.
+The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true
+(\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}).
+This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
+\[
+  z = s_{min} + C (s) (H - s_{min})
+  % \label{eq:SH94_1}
+\]
+where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and
+allows a $z$-coordinate to placed on top of the stretched coordinate,
+and $z$ is the depth (negative down from the asea surface).
+\begin{gather*}
+  s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1
+  % \label{eq:DOM_s}
+ \\
+  % \label{eq:DOM_sco_function}
+  C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)}
+\end{gather*}
+A stretching function,
+modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}),
+is also available and is more commonly used for shelf seas modelling:
+\[
+  C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)}
+         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)}
+                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)}
+  % \label{eq:SH94_2}
+\]
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]
+  \begin{center}
+    \includegraphics[]{Fig_sco_function}
+    \caption{
+      \protect\label{fig:sco_function}
+      Examples of the stretching function applied to a seamount;
+      from left to right: surface, surface and bottom, and bottom intensified resolutions
+    }
+  \end{center}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to
+the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and
+bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$.
+$b$ has been designed to allow surface and/or bottom increase of the vertical resolution
+(\autoref{fig:sco_function}).
+Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in
+an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}.
+In this case the a stretching function $\gamma$ is defined such that:
+\begin{equation}
+  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1
+  % \label{eq:z}
+\end{equation}
+The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
+\begin{gather*}
+  % \label{eq:DOM_gamma_deriv}
+  \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt)
+           + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\
+  \intertext{Where:}
+  % \label{eq:DOM_gamma}
+  f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2}
+  \quad \text{and} \quad \sigma = \frac{k}{n - 1}
+\end{gather*}
+This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
+the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards
+the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and
+user prescribed surface (\np{rn\_zs}) and bottom depths.
+The bottom cell depth in this example is given as a function of water depth:
+\[
+  % \label{eq:DOM_zb}
+  Z_b = h a + b
+\]
+where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]
+  \includegraphics[]{Fig_DOM_compare_coordinates_surface}
+  \caption{
+    A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines),
+    a 50 level $Z$-coordinate (contoured surfaces) and
+    the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for
+    a idealised bathymetry that goes from $50~m$ to $5500~m$ depth.
+    For clarity every third coordinate surface is shown.
+  }
+  \label{fig:fig_compare_coordinates_surface}
+\end{figure}
+ % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
+This gives a smooth analytical stretching in computational space that is constrained to
+given specified surface and bottom grid cell thicknesses in real space.
+This is not to be confused with the hybrid schemes that
+superimpose geopotential coordinates on terrain following coordinates thus
+creating a non-analytical vertical coordinate that
+therefore may suffer from large gradients in the vertical resolutions.
+This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching,
+but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.
+As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than
+the critical depth $h_c$.
+In this example two options are available in depths shallower than $h_c$,
+with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false
+(the z-coordinate being equal to the depths of the stretched coordinate at $h_c$).
+Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as
+large slopes lead to hydrostatic consistency.
+A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented,
+and is output as part of the model mesh file at the start of the run.
+% -------------------------------------------------------------------------------------------------------------
+%        z*- or s*-coordinate
+% -------------------------------------------------------------------------------------------------------------
+\subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})}
+\label{subsec:DOM_zgr_star}
+This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.
+%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
+Additionally, hybrid combinations of the three main coordinates are available:
+$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e).
+A further choice related to vertical coordinate concerns the presence (or not) of ocean
+cavities beneath ice shelves within the model domain.  A setting of \np{ln\_isfcav} as
+\forcode{.true.} indicates that the domain contains  ocean cavities, otherwise the top,
+wet layer of the ocean will always be at the ocean surface.  This option is currently only
+available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied
+at the ocean/ice shelf interface.
+Within the model, the arrays describing the grid point depths and vertical scale factors
+are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before},
+\textit{now} and \textit{after} time step.  The time at which they are defined is
+indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.  They are updated at each
+model time step. The initial fixed reference coordinate system is held in variable names
+with a $\_0$ suffix.  When the linear free surface option is used
+(\np{ln\_linssh}\forcode{ = .true.}), \textit{before}, \textit{now} and \textit{after}
+arrays are initially set to their reference counterpart and remain fixed.
+\subsubsection{Required fields}
+\label{sec:DOM_zgr_fields}
+The explicit specification of a range of fields related to the vertical grid are required for the definition of a configuration. These include:
+\begin{Verbatim}[fontsize=\tiny]
+int    ln_zco, ln_zps, ln_sco            /* flags for z-coord, z-coord with partial steps and s-coord    */
+int    ln_isfcav                         /* flag  for ice shelf cavities                                 */
+double e3t_1d, e3w_1d                    /* reference vertical scale factors at T and W points           */
+double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */
+double e3uw_0, e3vw_0                    /* vertical scale factors 3D coordinate at UW and VW points     */
+int    bottom_level, top_level           /* last wet T-points, 1st wet T-points (for ice shelf cavities) */
+                                         /* For reference:                                               */
+float  bathy_metry                       /* bathymetry used in setting top and bottom levels             */
+\end{Verbatim}
+This set of vertical metrics is sufficient to describe the initial depth and thickness of
+every gridcell in the model regardless of the choice of vertical coordinate. With constant
+z-levels, e3 metrics will be uniform across each horizontal level. In the partial step
+case each e3 at the \np{bottom\_level} (and, possibly, \np{top\_level} if ice cavities are
+present) may vary from its horizontal neighbours. And, in s-coordinates, variations can
+occur throughout the water column. With the non-linear free-surface, all the coordinates
+behave more like the s-coordinate in that variations occurr throughout the water column
+with displacements related to the sea surface height. These variations are typically much
+smaller than those arising from bottom fitted coordinates. The values for vertical metrics
+supplied in the domain configuration file can be considered as those arising from a flat
+sea surface with zero elevation.
+The \np{bottom\_level} and \np{top\_level} 2D arrays define the \np{bottom\_level} and top
+wet levels in each grid column. Without ice cavities, \np{top\_level} is essentially a land
+mask (0 on land; 1 everywhere else). With ice cavities, \np{top\_level} determines the
+first wet point below the overlying ice shelf.
 % -------------------------------------------------------------------------------------------------------------
 %        level bathymetry and mask
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Level bathymetry and mask}
+\subsubsection{Level bathymetry and mask}
 \label{subsec:DOM_msk}
+Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with
+steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}.
+The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which
+gives the number of ocean levels (\ie those that are not masked) at each $t$-point.
+mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which
+gdept $\leq$ bathy.
+Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after
+mbathy is computed.
+Isolated grid points that do not communicate with another ocean point at the same level are eliminated.
+As for the representation of bathymetry, a 2D integer array, misfdep, is created.
+misfdep defines the level of the first wet $t$-point.
+All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked.
+By default, $misfdep(:,:) = 1$ and no cells are masked.
+In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into
+the cavities are performed in the \textit{zgr\_isf} routine.
+The compatibility between ice shelf draft and bathymetry is checked.
+All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked).
+If only one cell on the water column is opened at $t$-, $u$- or $v$-points,
+the bathymetry or the ice shelf draft is dug to fit this constrain.
+If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.
+From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows:
+From \np{top\_level} and \np{bottom\_level} fields, the mask fields are defined as follows:
 \begin{alignat*}{2}
   tmask(i,j,k) &= &  &
     \begin{cases}
 &\text{if $                  k  <    misfdep(i,j)$} \\
 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\
 &\text{if $                  k  >     mbathy(i,j)$}
+&\text{if $                  k  <    top\_level(i,j)$} \\
+&\text{if $bottom\_level(i,j) \leq k \leq   top\_level(i,j)$} \\
+&\text{if $                  k  >     bottom\_level(i,j)$}
     \end{cases}
   \\
 …
 exactly in the same way as for the bottom boundary.
+The specification of closed lateral boundaries requires that at least
+the first and last rows and columns of the \textit{mbathy} array are set to zero.
+In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
+the second one and its first column equal to the last but one (and so too the mask arrays)
+(see \autoref{fig:LBC_jperio}).
+%% The specification of closed lateral boundaries requires that at least
+%% the first and last rows and columns of the \textit{mbathy} array are set to zero.
+%% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
+%% the second one and its first column equal to the last but one (and so too the mask arrays)
+%% (see \autoref{fig:LBC_jperio}).
+%-------------------------------------------------------------------------------------------------
+%        Closed seas
+%-------------------------------------------------------------------------------------------------
+\subsection{Closed seas} \label{subsec:DOM_closea}
+When a global ocean is coupled to an atmospheric model it is better to represent all large
+water bodies (\eg great lakes, Caspian sea...) even if the model resolution does not allow
+their communication with the rest of the ocean.  This is unnecessary when the ocean is
+forced by fixed atmospheric conditions, so these seas can be removed from the ocean
+domain.  The user has the option to set the bathymetry in closed seas to zero (see
+\autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater
+imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it
+should be noted here that a successful use of these options requires appropriate mask
+fields to be present in the domain configuration file. Among the possibilities are:
+\begin{Verbatim}[fontsize=\tiny]
+int    closea_mask          /* non-zero values in closed sea areas for optional masking                  */
+int    closea_mask_rnf      /* non-zero values in closed sea areas with runoff locations (precip only)   */
+int    closea_mask_emp      /* non-zero values in closed sea areas with runoff locations (total emp)     */
+\end{Verbatim}
+% -------------------------------------------------------------------------------------------------------------
+%        Grid files
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Output grid files}
+\label{subsec:DOM_meshmask}
+\nlst{namcfg}
+Most of the arrays relating to a particular ocean model configuration dicussed in this
+chapter (grid-point position, scale factors) can be saved in a file if namelist parameter
+\np{ln\_write\_cfg} (namelist \ngn{namcfg}) is set to \forcode{.true.}; the output
+filename is set thorugh parameter \np{cn\_domcfg\_out}. This is only really useful
+if the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and
+checking or confirmation is required.
+\nlst{namdom}
+Alternatively, all the arrays relating to a particular ocean model configuration
+(grid-point position, scale factors, depths and masks) can be saved in a file called
+\texttt{mesh\_mask} if namelist parameter \np{ln\_meshmask} (namelist \ngn{namdom}) is set
+to \forcode{.true.}. This file contains additional fields that can be useful for
+post-processing applications
 % ================================================================
 % Domain: Initial State (dtatsd & istate)
 % ================================================================
+\section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
+\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]
+{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
 \label{sec:DTA_tsd}
 %-----------------------------------------namtsd-------------------------------------------
 \nlst{namtsd}
 %------------------------------------------------------------------------------------------
+Options are defined in \ngn{namtsd}.
+By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and
+salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
+Basic initial state options are defined in \ngn{namtsd}.  By default, the ocean starts
+from rest (the velocity field is set to zero) and the initialization of temperature and
+salinity fields is controlled through the \np{ln\_tsd\_init} namelist parameter.
 \begin{description}
 \item[\np{ln\_tsd\_init}~\forcode{= .true.}]
   use a T and S input files that can be given on the model grid itself or on their native input data grid.
   In the latter case,
   the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid
   (see \autoref{subsec:SBC_iof}).
   The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.
+  The computation is done in the \mdl{dtatsd} module.
+\item[\np{ln\_tsd\_init}~\forcode{= .false.}]
   use constant salinity value of $35.5~psu$ and an analytical profile of temperature
   (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
+\item[\np{ln\_tsd\_init}\forcode{= .true.}]
+  Use T and S input files that can be given on the model grid itself or on their native
+  input data grids.  In the latter case, the data will be interpolated on-the-fly both in
+  the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}).  The
+  information relating to the input files are specified in the \np{sn\_tem} and
+  \np{sn\_sal} structures.  The computation is done in the \mdl{dtatsd} module.
+\item[\np{ln\_tsd\_init}\forcode{= .false.}]
+  Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine
+  contained in \mdl{userdef\_istate}. The default version sets horizontally uniform T and
+  profiles as used in the  GYRE configuration (see \autoref{sec:CFG_gyre}).
 \end{description}

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DYN.tex

-                      r10499
+                      r11422
 %           Horizontal divergence and relative vorticity
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
+\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]
+{Horizontal divergence and relative vorticity (\protect\mdl{divcur})}
 \label{subsec:DYN_divcur}
 …
 %           Sea Surface Height evolution
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
+\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]
+{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})}
 \label{subsec:DYN_sshwzv}
 …
 Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to
 the sea surface height equation otherwise tracer content will not be conserved
 \citep{Griffies_al_MWR01, Leclair_Madec_OM09}.
+\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}.
 The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom,
 …
 %        Vorticity term
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Vorticity term (\protect\mdl{dynvor})}
+\subsection[Vorticity term (\textit{dynvor.F90})]
+{Vorticity term (\protect\mdl{dynvor})}
 \label{subsec:DYN_vor}
 %------------------------------------------nam_dynvor----------------------------------------------------
 …
 %                 enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
+\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens = .true.})]
+{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})}
 \label{subsec:DYN_vor_ens}
 …
 %                 energy conserving scheme
 %-------------------------------------------------------------
+\subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
+\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene = .true.})]
+{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})}
 \label{subsec:DYN_vor_ene}
 …
 %                 mix energy/enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) }
+\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix = .true.})]
+{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.})}
 \label{subsec:DYN_vor_mix}
 …
 %                 energy and enstrophy conserving scheme
 %-------------------------------------------------------------
+\subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) }
+\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een = .true.})]
+{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.})}
 \label{subsec:DYN_vor_een}
 …
 Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....}
 A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}.
+A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}.
 The idea is to get rid of the double averaging by considering triad combinations of vorticity.
 It is noteworthy that this solution is conceptually quite similar to the one proposed by
 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
 The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified
 for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.
+\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}).
+The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified
+for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme.
 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:
 \[
 …
 \begin{figure}[!ht]
   \begin{center}
     \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad}
+    \includegraphics[width=\textwidth]{Fig_DYN_een_triad}
     \caption{
       \protect\label{fig:DYN_een_triad}
 …
 (with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry)
 that tends to reinforce the topostrophy of the flow
 (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.
+(\ie the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}.
 Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as
 …
 (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).
 Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of
 the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}.
+the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}.
 Furthermore, used in combination with a partial steps representation of bottom topography,
 it improves the interaction between current and topography,
 leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}.
+leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}.
 %--------------------------------------------------------------------------------------------------------------
 %           Kinetic Energy Gradient term
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})}
+\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]
+{Kinetic energy gradient term (\protect\mdl{dynkeg})}
 \label{subsec:DYN_keg}
 …
 %           Vertical advection term
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Vertical advection term (\protect\mdl{dynzad}) }
+\subsection[Vertical advection term (\textit{dynzad.F90})]
+{Vertical advection term (\protect\mdl{dynzad})}
 \label{subsec:DYN_zad}
 …
 When \np{ln\_dynzad\_zts}\forcode{ = .true.},
 a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term.
 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
+This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
 Note that in this case,
 a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability,
 …
 %           Coriolis plus curvature metric terms
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) }
+\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]
+{Coriolis plus curvature metric terms (\protect\mdl{dynvor})}
 \label{subsec:DYN_cor_flux}
 …
 %           Flux form Advection term
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Flux form advection term (\protect\mdl{dynadv}) }
+\subsection[Flux form advection term (\textit{dynadv.F90})]
+{Flux form advection term (\protect\mdl{dynadv})}
 \label{subsec:DYN_adv_flux}
 …
 a $2^{nd}$ order centered finite difference scheme, CEN2,
 or a $3^{rd}$ order upstream biased scheme, UBS.
 The latter is described in \citet{Shchepetkin_McWilliams_OM05}.
+The latter is described in \citet{shchepetkin.mcwilliams_OM05}.
 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.
 In flux form, the schemes differ by the choice of a space and time interpolation to define the value of
 …
 %                 2nd order centred scheme
 %-------------------------------------------------------------
+\subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}
+\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2 = .true.})]
+{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})}
 \label{subsec:DYN_adv_cen2}
 …
 %                 UBS scheme
 %-------------------------------------------------------------
+\subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}
+\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs = .true.})]
+{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})}
 \label{subsec:DYN_adv_ubs}
 …
 where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$.
 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error
 \citep{Shchepetkin_McWilliams_OM05}.
 The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 It is not a \emph{positive} scheme, meaning that false extrema are permitted.
 …
 while the second term, which is the diffusion part of the scheme,
 is evaluated using the \textit{before} velocity (forward in time).
 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme.
+This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme.
 Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by
 one coefficient.
 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme.
 …
 %           Hydrostatic pressure gradient term
 % ================================================================
+\section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
+\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]
+{Hydrostatic pressure gradient (\protect\mdl{dynhpg})}
 \label{sec:DYN_hpg}
 %------------------------------------------nam_dynhpg---------------------------------------------------
 …
 %           z-coordinate with full step
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}
+\subsection[Full step $Z$-coordinate (\forcode{ln_dynhpg_zco = .true.})]
+{Full step $Z$-coordinate (\protect\np{ln\_dynhpg\_zco}\forcode{ = .true.})}
 \label{subsec:DYN_hpg_zco}
 …
 %           z-coordinate with partial step
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}
+\subsection[Partial step $Z$-coordinate (\forcode{ln_dynhpg_zps = .true.})]
+{Partial step $Z$-coordinate (\protect\np{ln\_dynhpg\_zps}\forcode{ = .true.})}
 \label{subsec:DYN_hpg_zps}
 …
 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast number of papers
 (\eg, \citet{Song1998, Shchepetkin_McWilliams_OM05}).
+(\eg, \citet{song_MWR98, shchepetkin.mcwilliams_OM05}).
 A number of different pressure gradient options are coded but the ROMS-like,
 density Jacobian with cubic polynomial method is currently disabled whilst known bugs are under investigation.
 $\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
+$\bullet$ Traditional coding (see for example \citet{madec.delecluse.ea_JPO96}: (\np{ln\_dynhpg\_sco}\forcode{ = .true.})
 \begin{equation}
   \label{eq:dynhpg_sco}
 …
 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}\forcode{ = .true.})
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05}
+$\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{shchepetkin.mcwilliams_OM05}
 (\np{ln\_dynhpg\_djc}\forcode{ = .true.}) (currently disabled; under development)
 Note that expression \autoref{eq:dynhpg_sco} is commonly used when the variable volume formulation is activated
 (\key{vvl}) because in that case, even with a flat bottom,
 the coordinate surfaces are not horizontal but follow the free surface \citep{Levier2007}.
+the coordinate surfaces are not horizontal but follow the free surface \citep{levier.treguier.ea_rpt07}.
 The pressure jacobian scheme (\np{ln\_dynhpg\_prj}\forcode{ = .true.}) is available as
 an improved option to \np{ln\_dynhpg\_sco}\forcode{ = .true.} when \key{vvl} is active.
 …
 corresponds to the water replaced by the ice shelf.
 This top pressure is constant over time.
 A detailed description of this method is described in \citet{Losch2008}.\\
+A detailed description of this method is described in \citet{losch_JGR08}.\\
 The pressure gradient due to ocean load is computed using the expression \autoref{eq:dynhpg_sco} described in
 …
 %           Time-scheme
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .true./.false.})}
+\subsection[Time-scheme (\forcode{ln_dynhpg_imp = .{true,false}.})]
+{Time-scheme (\protect\np{ln\_dynhpg\_imp}\forcode{ = .\{true,false\}}.)}
 \label{subsec:DYN_hpg_imp}
 …
 the physical phenomenon that controls the time-step is internal gravity waves (IGWs).
 A semi-implicit scheme for doubling the stability limit associated with IGWs can be used
 \citep{Brown_Campana_MWR78, Maltrud1998}.
+\citep{brown.campana_MWR78, maltrud.smith.ea_JGR98}.
 It involves the evaluation of the hydrostatic pressure gradient as
 an average over the three time levels $t-\rdt$, $t$, and $t+\rdt$
 …
 % Surface Pressure Gradient
 % ================================================================
+\section{Surface pressure gradient (\protect\mdl{dynspg})}
+\section[Surface pressure gradient (\textit{dynspg.F90})]
+{Surface pressure gradient (\protect\mdl{dynspg})}
 \label{sec:DYN_spg}
 %-----------------------------------------nam_dynspg----------------------------------------------------
 …
 which imposes a very small time step when an explicit time stepping is used.
 Two methods are proposed to allow a longer time step for the three-dimensional equations:
 the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt}),
+the filtered free surface, which is a modification of the continuous equations (see \autoref{eq:PE_flt?}),
 and the split-explicit free surface described below.
 The extra term introduced in the filtered method is calculated implicitly,
 …
 % Explicit free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Explicit free surface (\protect\key{dynspg\_exp})}
+\subsection[Explicit free surface (\texttt{\textbf{key\_dynspg\_exp}})]
+{Explicit free surface (\protect\key{dynspg\_exp})}
 \label{subsec:DYN_spg_exp}
 …
 % Split-explict free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Split-explicit free surface (\protect\key{dynspg\_ts})}
+\subsection[Split-explicit free surface (\texttt{\textbf{key\_dynspg\_ts}})]
+{Split-explicit free surface (\protect\key{dynspg\_ts})}
 \label{subsec:DYN_spg_ts}
 %------------------------------------------namsplit-----------------------------------------------------------
 …
 The split-explicit free surface formulation used in \NEMO (\key{dynspg\_ts} defined),
 also called the time-splitting formulation, follows the one proposed by \citet{Shchepetkin_McWilliams_OM05}.
+also called the time-splitting formulation, follows the one proposed by \citet{shchepetkin.mcwilliams_OM05}.
 The general idea is to solve the free surface equation and the associated barotropic velocity equations with
 a smaller time step than $\rdt$, the time step used for the three dimensional prognostic variables
 …
 \begin{equation}
   \label{eq:BT_dyn}
   \frac{\partial {\rm \overline{{\bf U}}_h} }{\partial t}=
   -f\;{\rm {\bf k}}\times {\rm \overline{{\bf U}}_h}
   -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \rm {\overline{{\bf U}}_h} + \rm {\overline{\bf G}}
+  \frac{\partial {\mathrm \overline{{\mathbf U}}_h} }{\partial t}=
+  -f\;{\mathrm {\mathbf k}}\times {\mathrm \overline{{\mathbf U}}_h}
+  -g\nabla _h \eta -\frac{c_b^{\textbf U}}{H+\eta} \mathrm {\overline{{\mathbf U}}_h} + \mathrm {\overline{\mathbf G}}
 \end{equation}
 \[
   % \label{eq:BT_ssh}
   \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\rm{\bf \overline{U}}}_h \,} \right]+P-E
+  \frac{\partial \eta }{\partial t}=-\nabla \cdot \left[ {\left( {H+\eta } \right) \; {\mathrm{\mathbf \overline{U}}}_h \,} \right]+P-E
 \]
 % \end{subequations}
 where $\rm {\overline{\bf G}}$ is a forcing term held constant, containing coupling term between modes,
+where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing coupling term between modes,
 surface atmospheric forcing as well as slowly varying barotropic terms not explicitly computed to gain efficiency.
 The third term on the right hand side of \autoref{eq:BT_dyn} represents the bottom stress
 (see section \autoref{sec:ZDF_bfr}), explicitly accounted for at each barotropic iteration.
 Temporal discretization of the system above follows a three-time step Generalized Forward Backward algorithm
 detailed in \citet{Shchepetkin_McWilliams_OM05}.
+detailed in \citet{shchepetkin.mcwilliams_OM05}.
 AB3-AM4 coefficients used in \NEMO follow the second-order accurate,
 "multi-purpose" stability compromise as defined in \citet{Shchepetkin_McWilliams_Bk08}
+"multi-purpose" stability compromise as defined in \citet{shchepetkin.mcwilliams_ibk09}
 (see their figure 12, lower left).
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=0.7\textwidth]{Fig_DYN_dynspg_ts}
+    \includegraphics[width=\textwidth]{Fig_DYN_dynspg_ts}
     \caption{
       \protect\label{fig:DYN_dynspg_ts}
 …
 and time splitting not compatible.
 Advective barotropic velocities are obtained by using a secondary set of filtering weights,
 uniquely defined from the filter coefficients used for the time averaging (\citet{Shchepetkin_McWilliams_OM05}).
+uniquely defined from the filter coefficients used for the time averaging (\citet{shchepetkin.mcwilliams_OM05}).
 Consistency between the time averaged continuity equation and the time stepping of tracers is here the key to
 obtain exact conservation.
 …
 external gravity waves in idealized or weakly non-linear cases.
 Although the damping is lower than for the filtered free surface,
 it is still significant as shown by \citet{Levier2007} in the case of an analytical barotropic Kelvin wave.
+it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave.
 %>>>>>===============
 …
 the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
 We have tried various forms of such filtering,
 with the following method discussed in \cite{Griffies_al_MWR01} chosen due to
+with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
 its stability and reasonably good maintenance of tracer conservation properties (see ??).
 …
 % Filtered free surface formulation
 %--------------------------------------------------------------------------------------------------------------
+\subsection{Filtered free surface (\protect\key{dynspg\_flt})}
+\subsection[Filtered free surface (\texttt{\textbf{key\_dynspg\_flt}})]
+{Filtered free surface (\protect\key{dynspg\_flt})}
 \label{subsec:DYN_spg_fltp}
 The filtered formulation follows the \citet{Roullet_Madec_JGR00} implementation.
+The filtered formulation follows the \citet{roullet.madec_JGR00} implementation.
 The extra term introduced in the equations (see \autoref{subsec:PE_free_surface}) is solved implicitly.
 The elliptic solvers available in the code are documented in \autoref{chap:MISC}.
 …
   \[
     % \label{eq:spg_flt}
     \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
+    \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t}= {\mathrm {\mathbf M}}
     - g \nabla \left( \tilde{\rho} \ \eta \right)
     - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right)
 …
   where $T_c$, is a parameter with dimensions of time which characterizes the force,
   $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density,
   and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
+  and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient,
   non-linear and viscous terms in \autoref{eq:PE_dyn}.
 }   %end gmcomment
 …
 % Lateral diffusion term
 % ================================================================
+\section{Lateral diffusion term and operators (\protect\mdl{dynldf})}
+\section[Lateral diffusion term and operators (\textit{dynldf.F90})]
+{Lateral diffusion term and operators (\protect\mdl{dynldf})}
 \label{sec:DYN_ldf}
 %------------------------------------------nam_dynldf----------------------------------------------------
 …
 % ================================================================
 \subsection[Iso-level laplacian (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})]
             {Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})}
+\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap = .true.})]
+{Iso-level laplacian operator (\protect\np{ln\_dynldf\_lap}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_lap}
 …
   \left\{
     \begin{aligned}
       D_u^{l{\rm {\bf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
+      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
           \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
         {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
       D_v^{l{\rm {\bf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
+      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
           \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
         {A_f^{lm} \;e_{3f} \zeta } \right]
 …
 %           Rotated laplacian operator
 %--------------------------------------------------------------------------------------------------------------
 \subsection[Rotated laplacian (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})]
             {Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})}
+\subsection[Rotated laplacian (\forcode{ln_dynldf_iso = .true.})]
+{Rotated laplacian operator (\protect\np{ln\_dynldf\_iso}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_iso}
 …
 %           Iso-level bilaplacian operator
 %--------------------------------------------------------------------------------------------------------------
 \subsection[Iso-level bilaplacian (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})]
             {Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})}
+\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap = .true.})]
+{Iso-level bilaplacian operator (\protect\np{ln\_dynldf\_bilap}\forcode{ = .true.})}
 \label{subsec:DYN_ldf_bilap}
 …
 %           Vertical diffusion term
 % ================================================================
+\section{Vertical diffusion term (\protect\mdl{dynzdf})}
+\section[Vertical diffusion term (\textit{dynzdf.F90})]
+{Vertical diffusion term (\protect\mdl{dynzdf})}
 \label{sec:DYN_zdf}
 %----------------------------------------------namzdf------------------------------------------------------
 …
 There are two main options for wetting and drying code (wd):
 (a) an iterative limiter (il) and (b) a directional limiter (dl).
 The directional limiter is based on the scheme developed by \cite{WarnerEtal13} for RO
+The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO
 MS
 which was in turn based on ideas developed for POM by \cite{Oey06}. The iterative
+which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative
 limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
 and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
 …
 %   Iterative limiters
 %-----------------------------------------------------------------------------------------
 \subsection   [Directional limiter (\textit{wet\_dry})]
          {Directional limiter (\mdl{wet\_dry})}
+\subsection[Directional limiter (\textit{wet\_dry.F90})]
+{Directional limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_directional_limiter}
 The principal idea of the directional limiter is that
 …
 \cite{WarnerEtal13} state that in their scheme the velocity masks at the cell faces for the baroclinic
+\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic
 timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
 or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
 …
 %-----------------------------------------------------------------------------------------
 \subsection   [Iterative limiter (\textit{wet\_dry})]
          {Iterative limiter (\mdl{wet\_dry})}
+\subsection[Iterative limiter (\textit{wet\_dry.F90})]
+{Iterative limiter (\mdl{wet\_dry})}
 \label{subsec:DYN_wd_iterative_limiter}
 \subsubsection [Iterative flux limiter (\textit{wet\_dry})]
          {Iterative flux limiter (\mdl{wet\_dry})}
+\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]
+{Iterative flux limiter (\mdl{wet\_dry})}
 \label{subsubsec:DYN_wd_il_spg_limiter}
 …
 \end{equation}
 Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\it [Q: Why is
+Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is
 this necessary/desirable?]}. Substituting from (\ref{dyn_wd_continuity_coef}) gives an
 expression for the coefficient needed to multiply the outward flux at this cell in order
 …
 %      Surface pressure gradients
 %----------------------------------------------------------------------------------------
 \subsubsection   [Modification of surface pressure gradients (\textit{dynhpg})]
          {Modification of surface pressure gradients (\mdl{dynhpg})}
+\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]
+{Modification of surface pressure gradients (\mdl{dynhpg})}
 \label{subsubsec:DYN_wd_il_spg}
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht] \begin{center}
 \includegraphics[width=0.8\textwidth]{Fig_WAD_dynhpg}
+\includegraphics[width=\textwidth]{Fig_WAD_dynhpg}
 \caption{ \label{Fig_WAD_dynhpg}
 Illustrations of the three possible combinations of the logical variables controlling the
 …
 conditions.
 \subsubsection   [Additional considerations (\textit{usrdef\_zgr})]
          {Additional considerations (\mdl{usrdef\_zgr})}
+\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]
+{Additional considerations (\mdl{usrdef\_zgr})}
 \label{subsubsec:WAD_additional}
 …
 %      The WAD test cases
 %----------------------------------------------------------------------------------------
 \subsection   [The WAD test cases (\textit{usrdef\_zgr})]
          {The WAD test cases (\mdl{usrdef\_zgr})}
+\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]
+{The WAD test cases (\mdl{usrdef\_zgr})}
 \label{WAD_test_cases}
 …
 % Time evolution term
 % ================================================================
+\section{Time evolution term (\protect\mdl{dynnxt})}
+\section[Time evolution term (\textit{dynnxt.F90})]
+{Time evolution term (\protect\mdl{dynnxt})}
 \label{sec:DYN_nxt}

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LBC.tex

-                      r10614
+                      r11422
 % Boundary Condition at the Coast
 % ================================================================
+\section{Boundary condition at the coast (\protect\np{rn\_shlat})}
+\section[Boundary condition at the coast (\texttt{rn\_shlat})]
+{Boundary condition at the coast (\protect\np{rn\_shlat})}
 \label{sec:LBC_coast}
 %--------------------------------------------nam_lbc-------------------------------------------------------
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=0.90\textwidth]{Fig_LBC_uv}
+    \includegraphics[width=\textwidth]{Fig_LBC_uv}
     \caption{
       \protect\label{fig:LBC_uv}
 …
 \begin{figure}[!p]
   \begin{center}
     \includegraphics[width=0.90\textwidth]{Fig_LBC_shlat}
+    \includegraphics[width=\textwidth]{Fig_LBC_shlat}
     \caption{
       \protect\label{fig:LBC_shlat}
 …
 % Boundary Condition around the Model Domain
 % ================================================================
+\section{Model domain boundary condition (\protect\np{jperio})}
+\section[Model domain boundary condition (\texttt{jperio})]
+{Model domain boundary condition (\protect\np{jperio})}
 \label{sec:LBC_jperio}
 …
 %        Closed, cyclic (\np{jperio}\forcode{ = 0..2})
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Closed, cyclic (\protect\np{jperio}\forcode{= [0127]})}
+\subsection[Closed, cyclic (\forcode{jperio = [0127]})]
+{Closed, cyclic (\protect\np{jperio}\forcode{ = [0127]})}
 \label{subsec:LBC_jperio012}
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_LBC_jperio}
+    \includegraphics[width=\textwidth]{Fig_LBC_jperio}
     \caption{
       \protect\label{fig:LBC_jperio}
 …
 %        North fold (\textit{jperio = 3 }to $6)$
 % -------------------------------------------------------------------------------------------------------------
+\subsection{North-fold (\protect\np{jperio}\forcode{ = 3..6})}
+\subsection[North-fold (\forcode{jperio = [3-6]})]
+{North-fold (\protect\np{jperio}\forcode{ = [3-6]})}
 \label{subsec:LBC_north_fold}
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=0.90\textwidth]{Fig_North_Fold_T}
+    \includegraphics[width=\textwidth]{Fig_North_Fold_T}
     \caption{
       \protect\label{fig:North_Fold_T}
 …
 % Exchange with neighbouring processors
 % ====================================================================
+\section{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})}
+\section[Exchange with neighbouring processors (\textit{lbclnk.F90}, \textit{lib\_mpp.F90})]
+{Exchange with neighbouring processors (\protect\mdl{lbclnk}, \protect\mdl{lib\_mpp})}
 \label{sec:LBC_mpp}
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=0.90\textwidth]{Fig_mpp}
+    \includegraphics[width=\textwidth]{Fig_mpp}
     \caption{
       \protect\label{fig:mpp}
 …
 \begin{figure}[!ht]
   \begin{center}
     \includegraphics[width=0.90\textwidth]{Fig_mppini2}
+    \includegraphics[width=\textwidth]{Fig_mppini2}
     \caption {
       \protect\label{fig:mppini2}
 …
 The BDY module was modelled on the OBC module (see NEMO 3.4) and shares many features and
 a similar coding structure \citep{Chanut2005}.
+a similar coding structure \citep{chanut_rpt05}.
 The specification of the location of the open boundary is completely flexible and
 allows for example the open boundary to follow an isobath or other irregular contour.
 …
 \label{subsec:BDY_FRS_scheme}
 The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95},
+The Flow Relaxation Scheme (FRS) \citep{davies_QJRMS76,engedahl_T95},
 applies a simple relaxation of the model fields to externally-specified values over
 a zone next to the edge of the model domain.
 …
 \label{subsec:BDY_flather_scheme}
 The \citet{Flather_JPO94} scheme is a radiation condition on the normal,
+The \citet{flather_JPO94} scheme is a radiation condition on the normal,
 depth-mean transport across the open boundary.
 It takes the form
 …
 \label{subsec:BDY_orlanski_scheme}
 The Orlanski scheme is based on the algorithm described by \citep{Marchesiello2001}, hereafter MMS.
+The Orlanski scheme is based on the algorithm described by \citep{marchesiello.mcwilliams.ea_OM01}, hereafter MMS.
 The adaptive Orlanski condition solves a wave plus relaxation equation at the boundary:
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_LBC_bdy_geom}
+    \includegraphics[width=\textwidth]{Fig_LBC_bdy_geom}
     \caption {
       \protect\label{fig:LBC_bdy_geom}
 …
 These restrictions mean that data files used with versions of the
 model prior to Version 3.4 may not work with Version 3.4 onwards.
 A \fortran utility {\it bdy\_reorder} exists in the TOOLS directory which
+A \fortran utility {\itshape bdy\_reorder} exists in the TOOLS directory which
 will re-order the data in old BDY data files.
 …
 \begin{figure}[!t]
   \begin{center}
     \includegraphics[width=1.0\textwidth]{Fig_LBC_nc_header}
+    \includegraphics[width=\textwidth]{Fig_LBC_nc_header}
     \caption {
       \protect\label{fig:LBC_nc_header}

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex

-                      r10442
+                      r11422
 These three aspects of the lateral diffusion are set through namelist parameters
 (see the \ngn{nam\_traldf} and \ngn{nam\_dynldf} below).
 Note that this chapter describes the standard implementation of iso-neutral tracer mixing,
 and Griffies's implementation, which is used if \np{traldf\_grif}\forcode{ = .true.},
 is described in Appdx\autoref{apdx:triad}
 %-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
+Note that this chapter describes the standard implementation of iso-neutral tracer mixing.
+Griffies's implementation, which is used if \np{ln\_traldf\_triad}\forcode{ = .true.},
+is described in \autoref{apdx:triad}
+%-----------------------------------namtra_ldf - namdyn_ldf--------------------------------------------
 \nlst{namtra_ldf}
 …
 %--------------------------------------------------------------------------------------------------------------
+% ================================================================
+% Lateral Mixing Operator
+% ================================================================
+\section[Lateral mixing operators]
+{Lateral mixing operators}
+\label{sec:LDF_op}
+We remind here the different lateral mixing operators that can be used. Further details can be found in \autoref{subsec:TRA_ldf_op} and  \autoref{sec:DYN_ldf}.
+\subsection[No lateral mixing (\forcode{ln_traldf_OFF}, \forcode{ln_dynldf_OFF})]
+{No lateral mixing (\protect\np{ln\_traldf\_OFF}, \protect\np{ln\_dynldf\_OFF})}
+It is possible to run without explicit lateral diffusion on tracers (\protect\np{ln\_traldf\_OFF}\forcode{ = .true.}) and/or
+momentum (\protect\np{ln\_dynldf\_OFF}\forcode{ = .true.}). The latter option is even recommended if using the
+UBS advection scheme on momentum (\np{ln\_dynadv\_ubs}\forcode{ = .true.},
+see \autoref{subsec:DYN_adv_ubs}) and can be useful for testing purposes.
+\subsection[Laplacian mixing (\forcode{ln_traldf_lap}, \forcode{ln_dynldf_lap})]
+{Laplacian mixing (\protect\np{ln\_traldf\_lap}, \protect\np{ln\_dynldf\_lap})}
+Setting \protect\np{ln\_traldf\_lap}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_lap}\forcode{ = .true.} enables
+a second order diffusion on tracers and momentum respectively. Note that in \NEMO 4, one can not combine
+Laplacian and Bilaplacian operators for the same variable.
+\subsection[Bilaplacian mixing (\forcode{ln_traldf_blp}, \forcode{ln_dynldf_blp})]
+{Bilaplacian mixing (\protect\np{ln\_traldf\_blp}, \protect\np{ln\_dynldf\_blp})}
+Setting \protect\np{ln\_traldf\_blp}\forcode{ = .true.} and/or \protect\np{ln\_dynldf\_blp}\forcode{ = .true.} enables
+a fourth order diffusion on tracers and momentum respectively. It is implemented by calling the above Laplacian operator twice.
+We stress again that from \NEMO 4, the simultaneous use Laplacian and Bilaplacian operators is not allowed.
 % ================================================================
 % Direction of lateral Mixing
 % ================================================================
+\section{Direction of lateral mixing (\protect\mdl{ldfslp})}
+\section[Direction of lateral mixing (\textit{ldfslp.F90})]
+{Direction of lateral mixing (\protect\mdl{ldfslp})}
 \label{sec:LDF_slp}
 …
 \gmcomment{
   we should emphasize here that the implementation is a rather old one.
   Better work can be achieved by using \citet{Griffies_al_JPO98, Griffies_Bk04} iso-neutral scheme.
+  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme.
+}
 …
 %gm%  caution I'm not sure the simplification was a good idea!
 These slopes are computed once in \rou{ldfslp\_init} when \np{ln\_sco}\forcode{ = .true.}rue,
+These slopes are computed once in \rou{ldf\_slp\_init} when \np{ln\_sco}\forcode{ = .true.},
 and either \np{ln\_traldf\_hor}\forcode{ = .true.} or \np{ln\_dynldf\_hor}\forcode{ = .true.}.
 …
 %In practice, \autoref{eq:ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth.
 %By definition, neutral surfaces are tangent to the local $in situ$ density \citep{McDougall1987}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
+%By definition, neutral surfaces are tangent to the local $in situ$ density \citep{mcdougall_JPO87}, therefore in \autoref{eq:ldfslp_iso}, all the derivatives have to be evaluated at the same local pressure (which in decibars is approximated by the depth in meters).
 %In the $z$-coordinate, the derivative of the  \autoref{eq:ldfslp_iso} numerator is evaluated at the same depth \nocite{as what?} ($T$-level, which is the same as the $u$- and $v$-levels), so  the $in situ$ density can be used for its evaluation.
 …
   thus the $in situ$ density can be used.
   This is not the case for the vertical derivatives: $\delta_{k+1/2}[\rho]$ is replaced by $-\rho N^2/g$,
   where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{McDougall1987}
+  where $N^2$ is the local Brunt-Vais\"{a}l\"{a} frequency evaluated following \citet{mcdougall_JPO87}
   (see \autoref{subsec:TRA_bn2}).
 …
 \item[$s$- or hybrid $s$-$z$- coordinate: ]
   in the current release of \NEMO, iso-neutral mixing is only employed for $s$-coordinates if
   the Griffies scheme is used (\np{traldf\_grif}\forcode{ = .true.};
   see Appdx \autoref{apdx:triad}).
+  the Griffies scheme is used (\np{ln\_traldf\_triad}\forcode{ = .true.};
+  see \autoref{apdx:triad}).
   In other words, iso-neutral mixing will only be accurately represented with a linear equation of state
   (\np{nn\_eos}\forcode{ = 1..2}).
+  (\np{ln\_seos}\forcode{ = .true.}).
   In the case of a "true" equation of state, the evaluation of $i$ and $j$ derivatives in \autoref{eq:ldfslp_iso}
   will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes.
 …
   Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for
   the constraint on iso-neutral fluxes.
   Following \citet{Griffies_Bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
+  Following \citet{griffies_bk04}, instead of specifying directly that there is a zero neutral diffusive flux of
   locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between
   the neutral direction diffusive fluxes of potential temperature and salinity:
 …
 This implementation is a rather old one.
 It is similar to the one proposed by Cox [1987], except for the background horizontal diffusion.
 Indeed, the Cox implementation of isopycnal diffusion in GFDL-type models requires
+It is similar to the one proposed by \citet{cox_OM87}, except for the background horizontal diffusion.
+Indeed, the \citet{cox_OM87} implementation of isopycnal diffusion in GFDL-type models requires
 a minimum background horizontal diffusion for numerical stability reasons.
 To overcome this problem, several techniques have been proposed in which the numerical schemes of
 the ocean model are modified \citep{Weaver_Eby_JPO97, Griffies_al_JPO98}.
 Griffies's scheme is now available in \NEMO if \np{traldf\_grif\_iso} is set true; see Appdx \autoref{apdx:triad}.
 Here, another strategy is presented \citep{Lazar_PhD97}:
+the ocean model are modified \citep{weaver.eby_JPO97, griffies.gnanadesikan.ea_JPO98}.
+Griffies's scheme is now available in \NEMO if \np{ln\_traldf\_triad}=\forcode{= .true.}; see \autoref{apdx:triad}.
+Here, another strategy is presented \citep{lazar_phd97}:
 a local filtering of the iso-neutral slopes (made on 9 grid-points) prevents the development of
 grid point noise generated by the iso-neutral diffusion operator (\autoref{fig:LDF_ZDF1}).
 …
 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
 contrary to the \citet{Griffies_al_JPO98} operator which has that property.
+contrary to the \citet{griffies.gnanadesikan.ea_JPO98} operator which has that property.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
   \begin{center}
     \includegraphics[width=0.70\textwidth]{Fig_LDF_ZDF1}
+    \includegraphics[width=\textwidth]{Fig_LDF_ZDF1}
     \caption {
       \protect\label{fig:LDF_ZDF1}
 …
 % In addition and also for numerical stability reasons \citep{Cox1987, Griffies_Bk04},
+% In addition and also for numerical stability reasons \citep{cox_OM87, griffies_bk04},
 % the slopes are bounded by $1/100$ everywhere. This limit is decreasing linearly
 % to zero fom $70$ meters depth and the surface (the fact that the eddies "feel" the
 % surface motivates this flattening of isopycnals near the surface).
 For numerical stability reasons \citep{Cox1987, Griffies_Bk04}, the slopes must also be bounded by
 $1/100$ everywhere.
+For numerical stability reasons \citep{cox_OM87, griffies_bk04}, the slopes must also be bounded by
+the namelist scalar \np{rn\_slpmax} (usually $1/100$) everywhere.
 This constraint is applied in a piecewise linear fashion, increasing from zero at the surface to
 $1/100$ at $70$ metres and thereafter decreasing to zero at the bottom of the ocean
 (the fact that the eddies "feel" the surface motivates this flattening of isopycnals near the surface).
+\colorbox{yellow}{The way slopes are tapered has be checked. Not sure that this is still what is actually done.}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht]
   \begin{center}
     \includegraphics[width=0.70\textwidth]{Fig_eiv_slp}
+    \includegraphics[width=\textwidth]{Fig_eiv_slp}
     \caption{
       \protect\label{fig:eiv_slp}
 …
 (see \autoref{sec:LBC_coast}).
-% ================================================================
-% Lateral Mixing Operator
-% ================================================================
-\section{Lateral mixing operators (\protect\mdl{traldf}, \protect\mdl{traldf}) }
-\label{sec:LDF_op}
 % ================================================================
 % Lateral Mixing Coefficients
 % ================================================================
+\section{Lateral mixing coefficient (\protect\mdl{ldftra}, \protect\mdl{ldfdyn}) }
+\section[Lateral mixing coefficient (\forcode{nn_aht_ijk_t}, \forcode{nn_ahm_ijk_t})]
+{Lateral mixing coefficient (\protect\np{nn\_aht\_ijk\_t}, \protect\np{nn\_ahm\_ijk\_t})}
 \label{sec:LDF_coef}
+Introducing a space variation in the lateral eddy mixing coefficients changes the model core memory requirement,
+adding up to four extra three-dimensional arrays for the geopotential or isopycnal second order operator applied to
+momentum.
+Six CPP keys control the space variation of eddy coefficients: three for momentum and three for tracer.
+The three choices allow:
+a space variation in the three space directions (\key{traldf\_c3d},  \key{dynldf\_c3d}),
+in the horizontal plane (\key{traldf\_c2d},  \key{dynldf\_c2d}),
+or in the vertical only (\key{traldf\_c1d},  \key{dynldf\_c1d}).
+The default option is a constant value over the whole ocean on both momentum and tracers.
+The number of additional arrays that have to be defined and the gridpoint position at which
+they are defined depend on both the space variation chosen and the type of operator used.
+The resulting eddy viscosity and diffusivity coefficients can be a function of more than one variable.
+Changes in the computer code when switching from one option to another have been minimized by
+introducing the eddy coefficients as statement functions
+(include file \textit{ldftra\_substitute.h90} and \textit{ldfdyn\_substitute.h90}).
+The functions are replaced by their actual meaning during the preprocessing step (CPP).
+The specification of the space variation of the coefficient is made in \mdl{ldftra} and \mdl{ldfdyn},
+or more precisely in include files \textit{traldf\_cNd.h90} and \textit{dynldf\_cNd.h90}, with N=1, 2 or 3.
+The user can modify these include files as he/she wishes.
+The way the mixing coefficient are set in the reference version can be briefly described as follows:
+\subsubsection{Constant mixing coefficients (default option)}
+When none of the \key{dynldf\_...} and \key{traldf\_...} keys are defined,
+a constant value is used over the whole ocean for momentum and tracers,
+which is specified through the \np{rn\_ahm0} and \np{rn\_aht0} namelist parameters.
+\subsubsection{Vertically varying mixing coefficients (\protect\key{traldf\_c1d} and \key{dynldf\_c1d})}
+The 1D option is only available when using the $z$-coordinate with full step.
+Indeed in all the other types of vertical coordinate,
+the depth is a 3D function of (\textbf{i},\textbf{j},\textbf{k}) and therefore,
+introducing depth-dependent mixing coefficients will require 3D arrays.
+In the 1D option, a hyperbolic variation of the lateral mixing coefficient is introduced in which
+the surface value is \np{rn\_aht0} (\np{rn\_ahm0}), the bottom value is 1/4 of the surface value,
+and the transition takes place around z=300~m with a width of 300~m
+(\ie both the depth and the width of the inflection point are set to 300~m).
+This profile is hard coded in file \textit{traldf\_c1d.h90}, but can be easily modified by users.
+\subsubsection{Horizontally varying mixing coefficients (\protect\key{traldf\_c2d} and \protect\key{dynldf\_c2d})}
+By default the horizontal variation of the eddy coefficient depends on the local mesh size and
+The specification of the space variation of the coefficient is made in modules \mdl{ldftra} and \mdl{ldfdyn}.
+The way the mixing coefficients are set in the reference version can be described as follows:
+\subsection[Mixing coefficients read from file (\forcode{nn_aht_ijk_t = -20, -30}, \forcode{nn_ahm_ijk_t = -20,-30})]
+{ Mixing coefficients read from file (\protect\np{nn\_aht\_ijk\_t}\forcode{ = -20, -30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = -20, -30})}
+Mixing coefficients can be read from file if a particular geographical variation is needed. For example, in the ORCA2 global ocean model,
+the laplacian viscosity operator uses $A^l$~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
+decreases linearly to $A^l$~= 2.10$^3$ m$^2$/s at the equator \citep{madec.delecluse.ea_JPO96, delecluse.madec_icol99}.
+Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of ORCA2 and ORCA05.
+The provided fields can either be 2d (\np{nn\_aht\_ijk\_t}\forcode{ = -20}, \np{nn\_ahm\_ijk\_t}\forcode{ = -20}) or 3d (\np{nn\_aht\_ijk\_t}\forcode{ = -30},  \np{nn\_ahm\_ijk\_t}\forcode{ = -30}). They must be given at U, V points for tracers and T, F points for momentum (see \autoref{tab:LDF_files}).
+%-------------------------------------------------TABLE---------------------------------------------------
+\begin{table}[htb]
+  \begin{center}
+    \begin{tabular}{|l|l|l|l|}
+      \hline
+      Namelist parameter                        & Input filename                               & dimensions & variable names                      \\  \hline
+      \np{nn\_ahm\_ijk\_t}\forcode{ = -20}       & \forcode{eddy_viscosity_2D.nc }            &     $(i,j)$         & \forcode{ahmt_2d, ahmf_2d}  \\  \hline
+      \np{nn\_aht\_ijk\_t}\forcode{ = -20}           & \forcode{eddy_diffusivity_2D.nc }           &     $(i,j)$          & \forcode{ahtu_2d, ahtv_2d}    \\   \hline
+      \np{nn\_ahm\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_viscosity_3D.nc }            &     $(i,j,k)$          & \forcode{ahmt_3d, ahmf_3d}  \\  \hline
+      \np{nn\_aht\_ijk\_t}\forcode{ = -30}       & \forcode{eddy_diffusivity_3D.nc }           &     $(i,j,k)$         & \forcode{ahtu_3d, ahtv_3d}    \\   \hline
+    \end{tabular}
+    \caption{
+      \protect\label{tab:LDF_files}
+      Description of expected input files if mixing coefficients are read from NetCDF files.
+    }
+  \end{center}
+\end{table}
+%--------------------------------------------------------------------------------------------------------------
+\subsection[Constant mixing coefficients (\forcode{nn_aht_ijk_t = 0}, \forcode{nn_ahm_ijk_t = 0})]
+{ Constant mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 0}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 0})}
+If constant, mixing coefficients are set thanks to a velocity and a length scales ($U_{scl}$, $L_{scl}$) such that:
+\begin{equation}
+  \label{eq:constantah}
+  A_o^l = \left\{
+    \begin{aligned}
+      & \frac{1}{2} U_{scl} L_{scl}            & \text{for laplacian operator } \\
+      & \frac{1}{12} U_{scl} L_{scl}^3                    & \text{for bilaplacian operator }
+    \end{aligned}
+  \right.
+\end{equation}
+ $U_{scl}$ and $L_{scl}$ are given by the namelist parameters \np{rn\_Ud}, \np{rn\_Uv}, \np{rn\_Ld} and \np{rn\_Lv}.
+\subsection[Vertically varying mixing coefficients (\forcode{nn_aht_ijk_t = 10}, \forcode{nn_ahm_ijk_t = 10})]
+{Vertically varying mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 10}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 10})}
+In the vertically varying case, a hyperbolic variation of the lateral mixing coefficient is introduced in which
+the surface value is given by \autoref{eq:constantah}, the bottom value is 1/4 of the surface value,
+and the transition takes place around z=500~m with a width of 200~m.
+This profile is hard coded in module \mdl{ldfc1d\_c2d}, but can be easily modified by users.
+\subsection[Mesh size dependent mixing coefficients (\forcode{nn_aht_ijk_t = 20}, \forcode{nn_ahm_ijk_t = 20})]
+{Mesh size dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 20}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 20})}
+In that case, the horizontal variation of the eddy coefficient depends on the local mesh size and
 the type of operator used:
 \begin{equation}
 …
   A_l = \left\{
     \begin{aligned}
       & \frac{\max(e_1,e_2)}{e_{max}} A_o^l           & \text{for laplacian operator } \\
       & \frac{\max(e_1,e_2)^{3}}{e_{max}^{3}} A_o^l          & \text{for bilaplacian operator }
+      & \frac{1}{2} U_{scl}  \max(e_1,e_2)         & \text{for laplacian operator } \\
+      & \frac{1}{12} U_{scl}  \max(e_1,e_2)^{3}             & \text{for bilaplacian operator }
     \end{aligned}
   \right.
 \end{equation}
+where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain,
+and $A_o^l$ is the \np{rn\_ahm0} (momentum) or \np{rn\_aht0} (tracer) namelist parameter.
+where $U_{scl}$ is a user defined velocity scale (\np{rn\_Ud}, \np{rn\_Uv}).
 This variation is intended to reflect the lesser need for subgrid scale eddy mixing where
 the grid size is smaller in the domain.
 It was introduced in the context of the DYNAMO modelling project \citep{Willebrand_al_PO01}.
 Note that such a grid scale dependance of mixing coefficients significantly increase the range of stability of
 model configurations presenting large changes in grid pacing such as global ocean models.
+It was introduced in the context of the DYNAMO modelling project \citep{willebrand.barnier.ea_PO01}.
+Note that such a grid scale dependance of mixing coefficients significantly increases the range of stability of
+model configurations presenting large changes in grid spacing such as global ocean models.
 Indeed, in such a case, a constant mixing coefficient can lead to a blow up of the model due to
 large coefficient compare to the smallest grid size (see \autoref{sec:STP_forward_imp}),
 especially when using a bilaplacian operator.
+Other formulations can be introduced by the user for a given configuration.
+For example, in the ORCA2 global ocean model (see Configurations),
+the laplacian viscosity operator uses \np{rn\_ahm0}~= 4.10$^4$ m$^2$/s poleward of 20$^{\circ}$ north and south and
+decreases linearly to \np{rn\_aht0}~= 2.10$^3$ m$^2$/s at the equator \citep{Madec_al_JPO96, Delecluse_Madec_Bk00}.
+This modification can be found in routine \rou{ldf\_dyn\_c2d\_orca} defined in \mdl{ldfdyn\_c2d}.
+Similar modified horizontal variations can be found with the Antarctic or Arctic sub-domain options of
+ORCA2 and ORCA05 (see \&namcfg namelist).
+\subsubsection{Space varying mixing coefficients (\protect\key{traldf\_c3d} and \key{dynldf\_c3d})}
+The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases,
+\colorbox{yellow}{CASE \np{nn\_aht\_ijk\_t} = 21 to be added}
+\subsection[Mesh size and depth dependent mixing coefficients (\forcode{nn_aht_ijk_t = 30}, \forcode{nn_ahm_ijk_t = 30})]
+{Mesh size and depth dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 30}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 30})}
+The 3D space variation of the mixing coefficient is simply the combination of the 1D and 2D cases above,
 \ie a hyperbolic tangent variation with depth associated with a grid size dependence of
 the magnitude of the coefficient.
+\subsubsection{Space and time varying mixing coefficients}
+There is no default specification of space and time varying mixing coefficient.
+The only case available is specific to the ORCA2 and ORCA05 global ocean configurations.
+It provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and
+eddy induced velocity (ORCA05) that depends on the local growth rate of baroclinic instability.
+This specification is actually used when an ORCA key and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
+\subsection[Velocity dependent mixing coefficients (\forcode{nn_aht_ijk_t = 31}, \forcode{nn_ahm_ijk_t = 31})]
+{Flow dependent mixing coefficients (\protect\np{nn\_aht\_ijk\_t}\forcode{ = 31}, \protect\np{nn\_ahm\_ijk\_t}\forcode{ = 31})}
+In that case, the eddy coefficient is proportional to the local velocity magnitude so that the Reynolds number $Re =  \lvert U \rvert  e / A_l$ is constant (and here hardcoded to $12$):
+\colorbox{yellow}{JC comment: The Reynolds is effectively set to 12 in the code for both operators but shouldn't it be 2 for Laplacian ?}
+\begin{equation}
+  \label{eq:flowah}
+  A_l = \left\{
+    \begin{aligned}
+      & \frac{1}{12} \lvert U \rvert e          & \text{for laplacian operator } \\
+      & \frac{1}{12} \lvert U \rvert e^3             & \text{for bilaplacian operator }
+    \end{aligned}
+  \right.
+\end{equation}
+\subsection[Deformation rate dependent viscosities (\forcode{nn_ahm_ijk_t = 32})]
+{Deformation rate dependent viscosities (\protect\np{nn\_ahm\_ijk\_t}\forcode{ = 32})}
+This option refers to the \citep{smagorinsky_MW63} scheme which is here implemented for momentum only. Smagorinsky chose as a
+characteristic time scale $T_{smag}$ the deformation rate and for the lengthscale $L_{smag}$ the maximum wavenumber possible on the horizontal grid, e.g.:
+\begin{equation}
+  \label{eq:smag1}
+  \begin{split}
+    T_{smag}^{-1} & = \sqrt{\left( \partial_x u - \partial_y v\right)^2 + \left( \partial_y u + \partial_x v\right)^2  } \\
+    L_{smag} & = \frac{1}{\pi}\frac{e_1 e_2}{e_1 + e_2}
+  \end{split}
+\end{equation}
+Introducing a user defined constant $C$ (given in the namelist as \np{rn\_csmc}), one can deduce the mixing coefficients as follows:
+\begin{equation}
+  \label{eq:smag2}
+  A_{smag} = \left\{
+    \begin{aligned}
+      & C^2  T_{smag}^{-1}  L_{smag}^2       & \text{for laplacian operator } \\
+      & \frac{C^2}{8} T_{smag}^{-1} L_{smag}^4        & \text{for bilaplacian operator }
+    \end{aligned}
+  \right.
+\end{equation}
+For stability reasons, upper and lower limits are applied on the resulting coefficient (see \autoref{sec:STP_forward_imp}) so that:
+\begin{equation}
+  \label{eq:smag3}
+    \begin{aligned}
+      & C_{min} \frac{1}{2}   \lvert U \rvert  e    < A_{smag} < C_{max} \frac{e^2}{   8\rdt}                 & \text{for laplacian operator } \\
+      & C_{min} \frac{1}{12} \lvert U \rvert  e^3 < A_{smag} < C_{max} \frac{e^4}{64 \rdt}                 & \text{for bilaplacian operator }
+    \end{aligned}
+\end{equation}
+where $C_{min}$ and $C_{max}$ are adimensional namelist parameters given by \np{rn\_minfac} and \np{rn\_maxfac} respectively.
+\subsection{About space and time varying mixing coefficients}
 The following points are relevant when the eddy coefficient varies spatially:
 …
 (\autoref{sec:dynldf_properties}).
-(3) for isopycnal diffusion on momentum or tracers, an additional purely horizontal background diffusion with
-uniform coefficient can be added by setting a non zero value of \np{rn\_ahmb0} or \np{rn\_ahtb0},
-a background horizontal eddy viscosity or diffusivity coefficient
-(namelist parameters whose default values are $0$).
-However, the technique used to compute the isopycnal slopes is intended to get rid of such a background diffusion,
-since it introduces spurious diapycnal diffusion (see \autoref{sec:LDF_slp}).
-(4) when an eddy induced advection term is used (\key{traldf\_eiv}),
-$A^{eiv}$, the eddy induced coefficient has to be defined.
-Its space variations are controlled by the same CPP variable as for the eddy diffusivity coefficient
-(\ie \key{traldf\_cNd}).
-(5) the eddy coefficient associated with a biharmonic operator must be set to a \emph{negative} value.
-(6) it is possible to use both the laplacian and biharmonic operators concurrently.

Context Navigation

Legend:

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

NEMO/branches/2019/fix_vvl_ticket1791/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

NEMO/branches/2019/fix_vvl_ticket1791/doc

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/HTML_latex2html.sh

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/main

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_A.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_B.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_E.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/annex_iso.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_ASM.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_CONFIG.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIA.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DIU.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DOM.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_DYN.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LBC.tex

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_LDF.tex