Changeset 3294 for trunk/DOC/TexFiles/Chapters

trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex

Property svn:executable deleted

r2541	r3294
8	8	\vspace{-40pt}
9	9
10		\small{
	10	{\small
11	11	The ocean engine of NEMO (Nucleus for European Modelling of the Ocean) is a primitive
12	12	equation model adapted to regional and global ocean circulation problems. It is intended to

trunk/DOC/TexFiles/Chapters/Annex_A.tex

Property svn:executable deleted

-                      r2282
+                      r3294
 % Chain rule
 % ================================================================
 \section{Chain rule of $s-$coordinate}
+\section{The chain rule for $s-$coordinates}
 \label{Apdx_A_continuity}
 In order to establish the set of Primitive Equation in curvilinear $s-$coordinates
+In order to establish the set of Primitive Equation in curvilinear $s$-coordinates
 ($i.e.$ an orthogonal curvilinear coordinate in the horizontal and an Arbitrary Lagrangian
 Eulerian (ALE) coordinate in the vertical), we start from the set of equations established
 …
 % continuity equation
 % ================================================================
 \section{Continuity Equation in $s-$coordinate}
+\section{Continuity Equation in $s-$coordinates}
 \label{Apdx_A_continuity}
 …
 Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
 Introducing the dia-surface velocity component, $\omega $, defined as
 the velocity relative to the moving $s-$surfaces and normal to them:
+the volume flux across the moving $s$-surfaces per unit horizontal area:
 \begin{equation} \label{Apdx_A_w_s}
 \omega  = w - w_s - \sigma _1 \,u - \sigma _2 \,v    \\
 …
 This formulation of the pressure gradient is characterised by the appearance of a term depending on the
 the sea surface height only (last term on the right hand side of expression \eqref{Apdx_A_grad_p}).
+This term will be abusively named \textit{surface pressure gradient} whereas the first term will be named
+This term will be loosely termed \textit{surface pressure gradient}
+whereas the first term will be termed the
 \textit{hydrostatic pressure gradient} by analogy to the $z$-coordinate formulation.
 In fact, the the true surface pressure gradient is $1/\rho_o \nabla (\rho \eta)$, and
 …
 To sum up, in a curvilinear $s$-coordinate system, the vector invariant momentum equation
 solved by the model has the same mathematical expression as the one in a curvilinear
 $z-$coordinate, but the pressure gradient term :
+$z-$coordinate, except for the pressure gradient term :
 \begin{subequations} \label{Apdx_A_dyn_vect}
 \begin{multline} \label{Apdx_A_PE_dyn_vect_u}
 …
 \end{subequations}
 Both formulation share the same hydrostatic pressure balance expressed in terms of
 hydrostatic pressure and density anmalies, $p_h'$ and $d=( \frac{\rho}{\rho_o}-1 )$:
+hydrostatic pressure and density anomalies, $p_h'$ and $d=( \frac{\rho}{\rho_o}-1 )$:
 \begin{equation} \label{Apdx_A_dyn_zph}
 \frac{\partial p_h'}{\partial k} = - d \, g \, e_3
 …
 It is important to realize that the change in coordinate system has only concerned
 the position on the vertical. It has not affected (\textbf{i},\textbf{j},\textbf{k}), the
 orthogonal curvilinear set of unit vector. ($u$,$v$) are always horizontal velocities
+orthogonal curvilinear set of unit vectors. ($u$,$v$) are always horizontal velocities
 so that their evolution is driven by \emph{horizontal} forces, in particular
 the pressure gradient. By contrast, $\omega$ is not $w$, the third component of the velocity,
 but the dia-surface velocity component, $i.e.$ the velocity relative to the moving
 $s-$surfaces and normal to them.
+but the dia-surface velocity component, $i.e.$ the volume flux across the moving
+$s$-surfaces per unit horizontal area.

trunk/DOC/TexFiles/Chapters/Annex_B.tex

Property svn:executable deleted

-                      r2282
+                      r3294
 \label{Apdx_B_1}
 In the $z$-coordinate, the horizontal/vertical second order tracer diffusion operator
+\subsubsection*{In z-coordinates}
+In $z$-coordinates, the horizontal/vertical second order tracer diffusion operator
 is given by:
 \begin{eqnarray} \label{Apdx_B1}
  &D^T = \frac{1}{e_1 \, e_2}      \left[
   \left. \frac{\partial}{\partial i} \left(  \frac{e_2}{e_1}A^{lT} \;\left. \frac{\partial T}{\partial i} \right|_z   \right)   \right|_z      \right.
+ &D^T = \frac{1}{e_1 \, e_2}      \left[
+  \left. \frac{\partial}{\partial i} \left(  \frac{e_2}{e_1}A^{lT} \;\left. \frac{\partial T}{\partial i} \right|_z   \right)   \right|_z      \right.
                        \left.
 + \left. \frac{\partial}{\partial j} \left(  \frac{e_1}{e_2}A^{lT} \;\left. \frac{\partial T}{\partial j} \right|_z   \right)   \right|_z      \right]
++ \left. \frac{\partial}{\partial j} \left(  \frac{e_1}{e_2}A^{lT} \;\left. \frac{\partial T}{\partial j} \right|_z   \right)   \right|_z      \right]
 + \frac{\partial }{\partial z}\left( {A^{vT} \;\frac{\partial T}{\partial z}} \right)
 \end{eqnarray}
+In the $s$-coordinate, we defined the slopes of $s$-surfaces, $\sigma_1$ and
+$\sigma_2$ by \eqref{Apdx_A_s_slope} and the vertical/horizontal ratio of diffusion
+\subsubsection*{In generalized vertical coordinates}
+In $s$-coordinates, we defined the slopes of $s$-surfaces, $\sigma_1$ and
+$\sigma_2$ by \eqref{Apdx_A_s_slope} and the vertical/horizontal ratio of diffusion
 coefficient by $\epsilon = A^{vT} / A^{lT}$. The diffusion operator is given by:
 \begin{equation} \label{Apdx_B2}
 D^T = \left. \nabla \right|_s \cdot
+D^T = \left. \nabla \right|_s \cdot
            \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\
 \;\;\text{where} \;\Re =\left( {{\begin{array}{*{20}c}
 \hfill & 0 \hfill & {-\sigma _1 } \hfill \\
 \hfill & 1 \hfill & {-\sigma _2 } \hfill \\
  {-\sigma _1 } \hfill & {-\sigma _2 } \hfill & {\varepsilon +\sigma _1
+ {-\sigma _1 } \hfill & {-\sigma _2 } \hfill & {\varepsilon +\sigma _1
 ^2+\sigma _2 ^2} \hfill \\
 \end{array} }} \right)
 …
 or in expanded form:
 \begin{subequations}
 \begin{align*} {\begin{array}{*{20}l}
 D^T=& \frac{1}{e_1\,e_2\,e_3 }\;\left[ {\ \ \ \ e_2\,e_3\,A^{lT} \;\left.
+\begin{align*} {\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]
 \end{array} }
+&\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]
+\end{array} }
 \end{align*}
 \end{subequations}
 Equation \eqref{Apdx_B2} is obtained from \eqref{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 Appendix~\ref{Apdx_A}
 and use \eqref{Apdx_A_s_slope} and \eqref{Apdx_A_s_chain_rule}.
 Since no cross horizontal derivative $\partial _i \partial _j $ appears in
 \eqref{Apdx_B1}, the ($i$,$z$) and ($j$,$z$) planes are independent.
 The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$)
+Equation \eqref{Apdx_B2} is obtained from \eqref{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 Appendix~\ref{Apdx_A}
+and use \eqref{Apdx_A_s_slope} and \eqref{Apdx_A_s_chain_rule}.
+Since no cross horizontal derivative $\partial _i \partial _j $ appears in
+\eqref{Apdx_B1}, the ($i$,$z$) and ($j$,$z$) planes are independent.
+The derivation can then be demonstrated for the ($i$,$z$)~$\to$~($j$,$s$)
 transformation without any loss of generality:
 \begin{subequations}
 \begin{align*} {\begin{array}{*{20}l}
 D^T&=\frac{1}{e_1\,e_2} \left. {\frac{\partial }{\partial i}\left( {\frac{e_2}{e_1}A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_z } \right)} \right|_z
+\begin{subequations}
+\begin{align*} {\begin{array}{*{20}l}
+D^T&=\frac{1}{e_1\,e_2} \left. {\frac{\partial }{\partial i}\left( {\frac{e_2}{e_1}A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_z } \right)} \right|_z
                      +\frac{\partial }{\partial z}\left( {A^{vT}\;\frac{\partial T}{\partial z}} \right)     \\
+ \\
+%
+&=\frac{1}{e_1\,e_2 }\left[ {\left. {\;\frac{\partial }{\partial i}\left( {\frac{e_2}{e_1}A^{lT}\;\left( {\left. {\frac{\partial T}{\partial i}} \right|_s
+                                                    -\frac{e_1\,\sigma _1 }{e_3 }\frac{\partial T}{\partial s}} \right)} \right)} \right|_s } \right. \\
+& \qquad \qquad \left. { -\frac{e_1\,\sigma _1 }{e_3 }\frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\left( {\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{e_1 \,\sigma _1 }{e_3 }\frac{\partial T}{\partial s}} \right)} \right|_s } \right)\;} \right]
+\shoveright{ +\frac{1}{e_3 }\frac{\partial }{\partial s}\left[ {\frac{A^{vT}}{e_3 }\;\frac{\partial T}{\partial s}} \right]}  \qquad \qquad \qquad \\
+ \\
+%
+&=\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{e_2 }{e_1}A^{lT}\;\frac{\partial e_3 }{\partial i}} \right|_s \left. {\frac{\partial T}{\partial i}} \right|_s } \right. \\
+&  \qquad \qquad \quad \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 -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{ -e_1 \,\sigma _1 \frac{\partial }{\partial s}\left( {-\frac{e_2 \,\sigma _1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)\;\,\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\quad} \right] }\\
+\end{array} }     \\
+%
+ {\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:}
+%
+& =\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] }\\
 \\
+\allowdisplaybreaks
+&=\frac{1}{e_1\,e_2 }\left[ {\left. {\;\frac{\partial }{\partial i}\left( {\frac{e_2}{e_1}A^{lT}\;\left( {\left. {\frac{\partial T}{\partial i}} \right|_s
+                                                    -\frac{e_1\,\sigma _1 }{e_3 }\frac{\partial T}{\partial s}} \right)} \right)} \right|_s } \right. \\
+& \qquad \qquad \left. { -\frac{e_1\,\sigma _1 }{e_3 }\frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\left( {\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{e_1 \,\sigma _1 }{e_3 }\frac{\partial T}{\partial s}} \right)} \right|_s } \right)\;} \right]
+\shoveright{ +\frac{1}{e_3 }\frac{\partial }{\partial s}\left[ {\frac{A^{vT}}{e_3 }\;\frac{\partial T}{\partial s}} \right]}  \qquad \qquad \qquad \\
+\\
+\allowdisplaybreaks
+&=\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{e_2 }{e_1}A^{lT}\;\frac{\partial e_3 }{\partial i}} \right|_s \left. {\frac{\partial T}{\partial i}} \right|_s } \right. \\
+&  \qquad \qquad \quad \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 -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{ -e_1 \,\sigma _1 \frac{\partial }{\partial s}\left( {-\frac{e_2 \,\sigma _1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)\;\,\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\quad} \right] }\\
+\end{array} }     \\
+ {\begin{array}{*{20}l}
+%
+\allowdisplaybreaks
+\intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma _1 }{\partial s}$, 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 } \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] }\\
+\\
+&=\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]} \\
+%
+\allowdisplaybreaks
+&=\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]}
+\end{array} } \\
+{\begin{array}{*{20}l}
+%
 \intertext{using the same remark as just above, 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.\;\;\; \\
+& \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] }\\
+%
+\allowdisplaybreaks
+\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
+&= \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] }
+ \end{array} } \\
+{\begin{array}{*{20}l}
+%
+\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:}
+%
 & =\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]} \\
+%
+\allowdisplaybreaks
 \intertext{in other words, the horizontal Laplacian operator in the ($i$,$s$) plane takes the following form :}
 \end{array} }     \\
+%
 D^T = {\frac{1}{e_1\,e_2\,e_3}}
+& =\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:}
+\end{array} } \\
+%
+{\frac{1}{e_1\,e_2\,e_3}}
 \left( {{\begin{array}{*{30}c}
 {\left. {\frac{\partial \left( {e_2 e_3 \bullet } \right)}{\partial i}} \right|_s } \hfill \\
 …
 \left( {{\begin{array}{*{30}c}
  {1} \hfill & {-\sigma_1 } \hfill \\
  {-\sigma_1} \hfill & {\varepsilon_1^2} \hfill \\
 \end{array} }} \right)
 \cdot
+ {-\sigma_1} \hfill & {\varepsilon + \sigma_1^2} \hfill \\
+\end{array} }} \right)
+\cdot
 \left( {{\begin{array}{*{30}c}
 {\frac{1}{e_1 }\;\left. {\frac{\partial \bullet }{\partial i}} \right|_s } \hfill \\
 …
 \end{align*}
 \end{subequations}
+\addtocounter{equation}{-2}
 % ================================================================
 …
 \label{Apdx_B_2}
+The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in the ($i$,$j$,$k$)
+curvilinear coordinate system in which the equations of the ocean circulation model are
+\subsubsection*{In z-coordinates}
+The iso/diapycnal diffusive tensor $\textbf {A}_{\textbf I}$ expressed in 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}:
 \begin{equation*}
+\begin{equation} \label{Apdx_B3}
 \textbf {A}_{\textbf I} = \frac{A^{lT}}{\left( {1+a_1 ^2+a_2 ^2} \right)}
 \left[ {{\begin{array}{*{20}c}
 …
  {-a_1 } \hfill & {-a_2 } \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:
+\end{equation}
+where ($a_1$, $a_2$) are the isopycnal slopes in ($\textbf{i}$,
+$\textbf{j}$) directions, relative to geopotentials:
 \begin{equation*}
 a_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
 \qquad , \qquad
 a_2 =\frac{e_3 }{e_2 }\left( {\frac{\partial \rho }{\partial j}}
+a_2 =\frac{e_3 }{e_2 }\left( {\frac{\partial \rho }{\partial j}}
 \right)\left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
 \end{equation*}
 In practice, the isopycnal slopes are generally less than $10^{-2}$ in the ocean, so
+In practice, isopycnal slopes are generally less than $10^{-2}$ in the ocean, so
 $\textbf {A}_{\textbf I}$ can be simplified appreciably \citep{Cox1987}:
+\begin{equation*}
+{\textbf{A}_{\textbf{I}}} \approx A^{lT}
+\begin{subequations} \label{Apdx_B4}
+\begin{equation} \label{Apdx_B4a}
+{\textbf{A}_{\textbf{I}}} \approx A^{lT}\;\Re\;\text{where} \;\Re =
 \left[ {{\begin{array}{*{20}c}
 \hfill & 0 \hfill & {-a_1 } \hfill \\
 \hfill & 1 \hfill & {-a_2 } \hfill \\
  {-a_1 } \hfill & {-a_2 } \hfill & {\varepsilon +a_1 ^2+a_2 ^2} \hfill \\
+\end{array} }} \right]
+\end{equation*}
+The resulting isopycnal operator conserves the quantity and dissipates its square.
+The demonstration of the first property is trivial as \eqref{Apdx_B2} is the divergence
+\end{array} }} \right],
+\end{equation}
+and the iso/dianeutral diffusive operator in $z$-coordinates is then
+\begin{equation}\label{Apdx_B4b}
+ D^T = \left. \nabla \right|_z \cdot
+           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_z T  \right]. \\
+\end{equation}
+\end{subequations}
+Physically, the full tensor \eqref{Apdx_B3}
+represents strong isoneutral diffusion on a plane parallel to the isoneutral
+surface and weak dianeutral diffusion perpendicular to this plane.
+However, the approximate `weak-slope' tensor \eqref{Apdx_B4a} represents strong
+diffusion along the isoneutral surface, with weak
+\emph{vertical}  diffusion -- the principal axes of the tensor are no
+longer orthogonal. This simplification also decouples
+the ($i$,$z$) and ($j$,$z$) planes of the tensor. The weak-slope operator therefore takes the same
+form, \eqref{Apdx_B4}, as \eqref{Apdx_B2}, the diffusion operator for geopotential
+diffusion written in non-orthogonal $i,j,s$-coordinates. Written out
+explicitly,
+\begin{multline} \label{Apdx_B_ldfiso}
+ D^T=\frac{1}{e_1 e_2 }\left\{
+ {\;\frac{\partial }{\partial i}\left[ {A_h \left( {\frac{e_2}{e_1}\frac{\partial T}{\partial i}-a_1 \frac{e_2}{e_3}\frac{\partial T}{\partial k}} \right)} \right]}
+ {+\frac{\partial}{\partial j}\left[ {A_h \left( {\frac{e_1}{e_2}\frac{\partial T}{\partial j}-a_2 \frac{e_1}{e_3}\frac{\partial T}{\partial k}} \right)} \right]\;} \right\} \\
+\shoveright{+\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A_h \left( {-\frac{a_1 }{e_1 }\frac{\partial T}{\partial i}-\frac{a_2 }{e_2 }\frac{\partial T}{\partial j}+\frac{\left( {a_1 ^2+a_2 ^2+\varepsilon} \right)}{e_3 }\frac{\partial T}{\partial k}} \right)} \right]}. \\
+\end{multline}
+The isopycnal diffusion operator \eqref{Apdx_B4},
+\eqref{Apdx_B_ldfiso} conserves tracer quantity and dissipates its
+square. The demonstration of the first property is trivial as \eqref{Apdx_B4} is the divergence
 of fluxes. Let us demonstrate the second one:
 \begin{equation*}
 \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
           = -\iiint\limits_D \nabla T\;.\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
+\iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
+          = -\iiint\limits_D \nabla T\;.\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv,
 \end{equation*}
 since
 \begin{subequations}
 \begin{align*} {\begin{array}{*{20}l}
 \nabla T\;.\left( {{\rm {\bf A}}_{\rm {\bf 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(
 {\frac{\partial T}{\partial j}} \right)^2} \right. \\
+and since
+\begin{subequations}
+\begin{align*} {\begin{array}{*{20}l}
+\nabla T\;.\left( {{\rm {\bf A}}_{\rm {\bf 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(
+{\frac{\partial T}{\partial j}} \right)^2} \right. \\
 &\qquad \qquad \qquad
+{ \left. -\,{2a_2 \frac{\partial T}{\partial j}\frac{\partial T}{\partial k}+\left( {a_1 ^2+a_2 ^2} \right)\left( {\frac{\partial T}{\partial k}} \right)^2} \right]} \\
+&=A_h \left[ {\left( {\frac{\partial T}{\partial i}-a_1 \frac{\partial T}{\partial k}} \right)^2+\left( {\frac{\partial T}{\partial j}-a_2 \frac{\partial T}{\partial k}} \right)^2} \right]      \\
+{ \left. -\,{2a_2 \frac{\partial T}{\partial j}\frac{\partial T}{\partial k}+\left( {a_1 ^2+a_2 ^2+\varepsilon} \right)\left( {\frac{\partial T}{\partial k}} \right)^2} \right]} \\
+&=A_h \left[ {\left( {\frac{\partial T}{\partial i}-a_1 \frac{\partial
+          T}{\partial k}} \right)^2+\left( {\frac{\partial T}{\partial
+          j}-a_2 \frac{\partial T}{\partial k}} \right)^2}
+  +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\
 & \geq 0
 \end{array} }
+\end{array} }
 \end{align*}
 \end{subequations}
+the property becomes obvious.
+The resulting diffusion operator in $z$-coordinate has the following form :
+\begin{multline*} \label{Apdx_B_ldfiso}
+ D^T=\frac{1}{e_1 e_2 }\left\{
+ {\;\frac{\partial }{\partial i}\left[ {A_h \left( {\frac{e_2}{e_1}\frac{\partial T}{\partial i}-a_1 \frac{e_2}{e_3}\frac{\partial T}{\partial k}} \right)} \right]}
+ {+\frac{\partial}{\partial j}\left[ {A_h \left( {\frac{e_1}{e_2}\frac{\partial T}{\partial j}-a_2 \frac{e_1}{e_3}\frac{\partial T}{\partial k}} \right)} \right]\;} \right\} \\
+\shoveright{+\frac{1}{e_3 }\frac{\partial }{\partial k}\left[ {A_h \left( {-\frac{a_1 }{e_1 }\frac{\partial T}{\partial i}-\frac{a_2 }{e_2 }\frac{\partial T}{\partial j}+\frac{\left( {a_1 ^2+a_2 ^2} \right)}{e_3 }\frac{\partial T}{\partial k}} \right)} \right]} \\
+\end{multline*}
+It has to be emphasised that the simplification introduced, leads to a decoupling
+between ($i$,$z$) and ($j$,$z$) planes. The operator has therefore the same
+expression as \eqref{Apdx_B3}, the diffusion operator obtained for geopotential
+diffusion in the $s$-coordinate.
+\addtocounter{equation}{-1}
+ the property becomes obvious.
+\subsubsection*{In generalized vertical coordinates}
+Because the weak-slope operator \eqref{Apdx_B4}, \eqref{Apdx_B_ldfiso} is decoupled
+in the ($i$,$z$) and ($j$,$z$) planes, it may be transformed into
+generalized $s$-coordinates in the same way as \eqref{Apdx_B_1} was transformed into
+\eqref{Apdx_B_2}. The resulting operator then takes the simple form
+\begin{equation} \label{Apdx_B_ldfiso_s}
+D^T = \left. \nabla \right|_s \cdot
+           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_s T  \right] \\
+\;\;\text{where} \;\Re =\left( {{\begin{array}{*{20}c}
+\hfill & 0 \hfill & {-r _1 } \hfill \\
+\hfill & 1 \hfill & {-r _2 } \hfill \\
+ {-r _1 } \hfill & {-r _2 } \hfill & {\varepsilon +r _1
+^2+r _2 ^2} \hfill \\
+\end{array} }} \right),
+\end{equation}
+where ($r_1$, $r_2$) are the isopycnal slopes in ($\textbf{i}$,
+$\textbf{j}$) directions, relative to $s$-coordinate surfaces:
+\begin{equation*}
+r_1 =\frac{e_3 }{e_1 }\left( {\frac{\partial \rho }{\partial i}} \right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1}
+\qquad , \qquad
+r_2 =\frac{e_3 }{e_2 }\left( {\frac{\partial \rho }{\partial j}}
+\right)\left( {\frac{\partial \rho }{\partial s}} \right)^{-1}.
+\end{equation*}
+To prove  \eqref{Apdx_B5}  by direct re-expression of \eqref{Apdx_B_ldfiso} is
+straightforward, but laborious. An easier way is first to note (by reversing the
+derivation of \eqref{Apdx_B_2} from \eqref{Apdx_B_1} ) that the
+weak-slope operator may be \emph{exactly} reexpressed in
+non-orthogonal $i,j,\rho$-coordinates as
+\begin{equation} \label{Apdx_B5}
+D^T = \left. \nabla \right|_\rho \cdot
+           \left[ A^{lT} \;\Re \cdot \left. \nabla \right|_\rho T  \right] \\
+\;\;\text{where} \;\Re =\left( {{\begin{array}{*{20}c}
+\hfill & 0 \hfill &0 \hfill \\
+\hfill & 1 \hfill & 0 \hfill \\
+\hfill & 0 \hfill & \varepsilon \hfill \\
+\end{array} }} \right).
+\end{equation}
+Then direct transformation from $i,j,\rho$-coordinates to
+$i,j,s$-coordinates gives \eqref{Apdx_B_ldfiso_s} immediately.
+Note that the weak-slope approximation is only made in
+transforming from 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,
+\eqref{Apdx_B_1} onto $s$-coordinates is exact, however steep the $s$-surfaces.
 % ================================================================
 …
 \label{Apdx_B_3}
 The second order momentum diffusion operator (Laplacian) in the $z$-coordinate
 is found by applying \eqref{Eq_PE_lap_vector}, the expression for the Laplacian
 of a vector,  to the horizontal velocity vector :
+The second order momentum diffusion operator (Laplacian) in the $z$-coordinate
+is found by applying \eqref{Eq_PE_lap_vector}, the expression for the Laplacian
+of a vector,  to the horizontal velocity vector :
 \begin{align*}
 \Delta {\textbf{U}}_h
+\Delta {\textbf{U}}_h
 &=\nabla \left( {\nabla \cdot {\textbf{U}}_h } \right)-
 \nabla \times \left( {\nabla \times {\textbf{U}}_h } \right)    \\
 …
  {\frac{1}{e_3 }\frac{\partial \chi }{\partial k}} \hfill \\
 \end{array} }} \right)-\left( {{\begin{array}{*{20}c}
  {\frac{1}{e_2 }\frac{\partial \zeta }{\partial j}-\frac{1}{e_3
 }\frac{\partial }{\partial k}\left( {\frac{1}{e_3 }\frac{\partial
+ {\frac{1}{e_2 }\frac{\partial \zeta }{\partial j}-\frac{1}{e_3
+}\frac{\partial }{\partial k}\left( {\frac{1}{e_3 }\frac{\partial
 u}{\partial k}} \right)} \hfill \\
  {\frac{1}{e_3 }\frac{\partial }{\partial k}\left( {-\frac{1}{e_3
 }\frac{\partial v}{\partial k}} \right)-\frac{1}{e_1 }\frac{\partial \zeta
+ {\frac{1}{e_3 }\frac{\partial }{\partial k}\left( {-\frac{1}{e_3
+}\frac{\partial v}{\partial k}} \right)-\frac{1}{e_1 }\frac{\partial \zeta
 }{\partial i}} \hfill \\
  {\frac{1}{e_1 e_2 }\left[ {\frac{\partial }{\partial i}\left( {\frac{e_2
 }{e_3 }\frac{\partial u}{\partial k}} \right)-\frac{\partial }{\partial
 j}\left( {-\frac{e_1 }{e_3 }\frac{\partial v}{\partial k}} \right)} \right]}
+ {\frac{1}{e_1 e_2 }\left[ {\frac{\partial }{\partial i}\left( {\frac{e_2
+}{e_3 }\frac{\partial u}{\partial k}} \right)-\frac{\partial }{\partial
+j}\left( {-\frac{e_1 }{e_3 }\frac{\partial v}{\partial k}} \right)} \right]}
 \hfill \\
 \end{array} }} \right)
 …
 \end{array} }} \right)
 \end{align*}
 Using \eqref{Eq_PE_div}, the definition of the horizontal divergence, the third
+Using \eqref{Eq_PE_div}, the definition of the horizontal divergence, the third
 componant of the second vector is obviously zero and thus :
 \begin{equation*}
 …
 \end{equation*}
 Note that this operator ensures a full separation between the vorticity and horizontal
 divergence fields (see Appendix~\ref{Apdx_C}). It is only equal to a Laplacian
+Note that this operator ensures a full separation between the vorticity and horizontal
+divergence fields (see Appendix~\ref{Apdx_C}). It is only equal to a Laplacian
 applied to each component in Cartesian coordinates, not on the sphere.
 The horizontal/vertical second order (Laplacian type) operator used to diffuse
+The horizontal/vertical second order (Laplacian type) operator used to diffuse
 horizontal momentum in the $z$-coordinate therefore takes the following form :
 \begin{equation} \label{Apdx_B_Lap_U}
  {\textbf{D}}^{\textbf{U}} =
+ {\textbf{D}}^{\textbf{U}} =
      \nabla _h \left( {A^{lm}\;\chi } \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 {\rm {\bf U}}_h }{\partial k}} \right) \\
 \end{equation}
 that is, in expanded form:
 \begin{align*}
 D^{\textbf{U}}_u
+D^{\textbf{U}}_u
 & = \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}      \\
 D^{\textbf{U}}_v
+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}
 …
 \end{align*}
 Note Bene: introducing a rotation in \eqref{Apdx_B_Lap_U} does not lead to a
 useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.
 Similarly, we did not found an expression of practical use for the geopotential
 horizontal/vertical Laplacian operator in the $s$-coordinate. Generally,
 \eqref{Apdx_B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is
+Note Bene: introducing a rotation in \eqref{Apdx_B_Lap_U} does not lead to a
+useful expression for the iso/diapycnal Laplacian operator in the $z$-coordinate.
+Similarly, we did not found an expression of practical use for the geopotential
+horizontal/vertical Laplacian operator in the $s$-coordinate. Generally,
+\eqref{Apdx_B_Lap_U} is used in both $z$- and $s$-coordinate systems, that is
 a Laplacian diffusion is applied on momentum along the coordinate directions.

trunk/DOC/TexFiles/Chapters/Annex_C.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Annex_D.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Annex_E.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Annex_ISO.tex

-                      r2376
+                      r3294
 % ================================================================
 % Iso-neutral diffusion :
+% Iso-neutral diffusion :
 % ================================================================
+\chapter{Griffies's iso-neutral diffusion}
+\label{Apdx_C}
+\chapter[Iso-neutral diffusion and eddy advection using
+triads]{Iso-neutral diffusion and eddy advection using triads}
+\label{sec:triad}
 \minitoc
+\section{Griffies's formulation of iso-neutral diffusion}
+\subsection{Introduction}
+ We define a scheme that get its inspiration from the scheme of
+\citet{Griffies_al_JPO98}, but is formulated within the \NEMO
+\pagebreak
+\section{Choice of namelist parameters}
+%-----------------------------------------nam_traldf------------------------------------------------------
+\namdisplay{namtra_ldf}
+%---------------------------------------------------------------------------------------------------------
+If the namelist variable \np{ln\_traldf\_grif} is set true (and
+\key{ldfslp} is set), \NEMO updates both active and passive tracers
+using the Griffies triad representation of iso-neutral diffusion and
+the eddy-induced advective skew (GM) fluxes. Otherwise (by default) the
+filtered version of Cox's original scheme is employed
+(\S\ref{LDF_slp}). In the present implementation of the Griffies
+scheme, the advective skew fluxes are implemented even if
+\key{traldf\_eiv} is not set.
+Values of iso-neutral diffusivity and GM coefficient are set as
+described in \S\ref{LDF_coef}. If none of the keys \key{traldf\_cNd},
+N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and
+GM diffusivity $A_e$ are directly set by \np{rn\_aeih\_0} and
+\np{rn\_aeiv\_0}. If 2D-varying coefficients are set with
+\key{traldf\_c2d} then $A_l$ is reduced in proportion with horizontal
+scale factor according to \eqref{Eq_title} \footnote{Except in global
+  $0.5^{\circ}$ runs (\key{orca\_r05}) with \key{traldf\_eiv}, where
+  $A_l$ is set like $A_e$ but with a minimum vale of
+  $100\;\mathrm{m}^2\;\mathrm{s}^{-1}$}. In idealised setups with
+\key{traldf\_c2d}, $A_e$ is reduced similarly, but if \key{traldf\_eiv}
+is set in the global configurations \key{orca\_r2}, \key{orca\_r1} or
+\key{orca\_r05} with \key{traldf\_c2d}, a horizontally varying $A_e$ is
+instead set from the Held-Larichev parameterisation\footnote{In this
+  case, $A_e$ at low latitudes $|\theta|<20^{\circ}$ is further
+  reduced by a factor $|f/f_{20}|$, where $f_{20}$ is the value of $f$
+  at $20^{\circ}$~N} (\mdl{ldfeiv}) and \np{rn\_aeiv\_0} is ignored
+unless it is zero.
+The options specific to the Griffies scheme include:
+\begin{description}[font=\normalfont]
+\item[\np{ln\_traldf\_gdia}] Default value is false. See \S\ref{sec:triad:sfdiag}. If this is set true, time-mean
+  eddy-advective (GM) velocities are output for diagnostic purposes, even
+  though the eddy advection is accomplished by means of the skew
+  fluxes.
+\item[\np{ln\_traldf\_iso}] See \S\ref{sec:triad:taper}. If this is set false (the default), then
+  `iso-neutral' mixing is accomplished within the surface mixed-layer
+  along slopes linearly decreasing with depth from the value immediately below
+  the mixed-layer to zero (flat) at the surface (\S\ref{sec:triad:lintaper}). This is the same
+  treatment as used in the default implementation
+  \S\ref{LDF_slp_iso}; Fig.~\ref{Fig_eiv_slp}.  Where
+  \np{ln\_traldf\_iso} is set true, 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 \S\ref{sec:triad:Gerdes-taper}
+\item[\np{ln\_traldf\_botmix}] See \S\ref{sec:triad:iso_bdry}. If this
+  is set false (the default) then the lateral diffusive fluxes
+  associated with triads partly masked by topography are neglected. If
+  it is set true, however, then these lateral diffusive fluxes are
+  applied, giving smoother bottom tracer fields at the cost of
+  introducing diapycnal mixing.
+\end{description}
+\section{Triad formulation of iso-neutral diffusion}
+\label{sec:triad:iso}
+We have implemented into \NEMO a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO
 framework, using scale factors rather than grid-sizes.
+The off-diagonal terms of the small angle diffusion tensor
+\eqref{Eq_PE_iso_tensor} produce skew-fluxes along the
+i- and j-directions resulting from the vertical tracer gradient:
+\begin{align}
+  \label{eq:i13c}
+  &+\kappa r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad+\kappa r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\
+\intertext{and in the k-direction resulting from the lateral tracer gradients}
+  \label{eq:i31c}
+ & \kappa r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\kappa r_2\frac{1}{e_1}\frac{\partial T}{\partial i}
+\end{align}
+where \eqref{Eq_PE_iso_slopes}
+\subsection{The iso-neutral diffusion operator}
+The iso-neutral second order tracer diffusive operator for small
+angles between iso-neutral surfaces and geopotentials is given by
+\eqref{Eq_PE_iso_tensor}:
+\begin{subequations} \label{eq:triad:PE_iso_tensor}
+  \begin{equation}
+    D^{lT}=-\Div\vect{f}^{lT}\equiv
+    -\frac{1}{e_1e_2e_3}\left[\pd{i}\left (f_1^{lT}e_2e_3\right) +
+      \pd{j}\left (f_2^{lT}e_2e_3\right) + \pd{k}\left (f_3^{lT}e_1e_2\right)\right],
+  \end{equation}
+  where the diffusive flux per unit area of physical space
+  \begin{equation}
+    \vect{f}^{lT}=-\Alt\Re\cdot\grad T,
+  \end{equation}
+  \begin{equation}
+    \label{eq:triad:PE_iso_tensor:c}
+    \mbox{with}\quad \;\;\Re =
+    \begin{pmatrix}
+&0&-r_1\mystrut \\
+&1&-r_2\mystrut \\
+      -r_1&-r_2&r_1 ^2+r_2 ^2\mystrut
+    \end{pmatrix}
+    \quad \text{and} \quad\grad T=
+    \begin{pmatrix}
+      \frac{1}{e_1}\pd[T]{i}\mystrut \\
+      \frac{1}{e_2}\pd[T]{j}\mystrut \\
+      \frac{1}{e_3}\pd[T]{k}\mystrut
+    \end{pmatrix}.
+  \end{equation}
+\end{subequations}
+% \left( {{\begin{array}{*{20}c}
+%  1 \hfill & 0 \hfill & {-r_1 } \hfill \\
+%  0 \hfill & 1 \hfill & {-r_2 } \hfill \\
+%  {-r_1 } \hfill & {-r_2 } \hfill & {r_1 ^2+r_2 ^2} \hfill \\
+% \end{array} }} \right)
+ Here \eqref{Eq_PE_iso_slopes}
 \begin{align*}
   r_1 &=-\frac{e_3 }{e_1 } \left( \frac{\partial \rho }{\partial i}
 …
     }{\partial k} \right)^{-1}
 \end{align*}
+is the i-component of the slope of the isoneutral surface relative to the computational
+surface, and $r_2$ is the j-component.
+The extra vertical diffusive flux associated with the $_{33}$
+is the $i$-component of the slope of the iso-neutral surface relative to the computational
+surface, and $r_2$ is the $j$-component.
+We will find it useful to consider the fluxes per unit area in $i,j,k$
+space; we write
+\begin{equation}
+  \label{eq:triad:Fijk}
+  \vect{F}_{\mathrm{iso}}=\left(f_1^{lT}e_2e_3, f_2^{lT}e_1e_3, f_3^{lT}e_1e_2\right).
+\end{equation}
+Additionally, we will sometimes write the contributions towards the
+fluxes $\vect{f}$ and $\vect{F}_\mathrm{iso}$ from the component
+$R_{ij}$ of $\Re$ as $f_{ij}$, $F_{\mathrm{iso}\: ij}$, with
+$f_{ij}=R_{ij}e_i^{-1}\partial T/\partial x_i$ (no summation) etc.
+The off-diagonal terms of the small angle diffusion tensor
+\eqref{Eq_PE_iso_tensor}, \eqref{eq:triad:PE_iso_tensor:c} produce skew-fluxes along the
+$i$- and $j$-directions resulting from the vertical tracer gradient:
+\begin{align}
+  \label{eq:triad:i13c}
+  f_{13}=&+\Alt r_1\frac{1}{e_3}\frac{\partial T}{\partial k},\qquad f_{23}=+\Alt r_2\frac{1}{e_3}\frac{\partial T}{\partial k}\\
+\intertext{and in the k-direction resulting from the lateral tracer gradients}
+  \label{eq:triad:i31c}
+ f_{31}+f_{32}=& \Alt r_1\frac{1}{e_1}\frac{\partial T}{\partial i}+\Alt r_2\frac{1}{e_1}\frac{\partial T}{\partial i}
+\end{align}
+The vertical diffusive flux associated with the $_{33}$
 component of the small angle diffusion tensor is
 \begin{equation}
   \label{eq:i33c}
   -\kappa(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.
+  \label{eq:triad:i33c}
+  f_{33}=-\Alt(r_1^2 +r_2^2) \frac{1}{e_3}\frac{\partial T}{\partial k}.
 \end{equation}
 Since there are no cross terms involving $r_1$ and $r_2$ in the above, we can
 consider the isoneutral diffusive fluxes separately in the i-k and j-k
+consider the iso-neutral diffusive fluxes separately in the $i$-$k$ and $j$-$k$
 planes, just adding together the vertical components from each
 plane. The following description will describe the fluxes on the i-k
+plane. The following description will describe the fluxes on the $i$-$k$
 plane.
 There is no natural discretization for the i-component of the
 skew-flux, \eqref{eq:i13c}, as
 although it must be evaluated at u-points, it involves vertical
+There is no natural discretization for the $i$-component of the
+skew-flux, \eqref{eq:triad:i13c}, as
+although it must be evaluated at $u$-points, it involves vertical
 gradients (both for the tracer and the slope $r_1$), defined at
 w-points. Similarly, the vertical skew flux, \eqref{eq:i31c}, is evaluated at
 w-points but involves horizontal gradients defined at u-points.
+$w$-points. Similarly, the vertical skew flux, \eqref{eq:triad:i31c}, is evaluated at
+$w$-points but involves horizontal gradients defined at $u$-points.
 \subsection{The standard discretization}
 The straightforward approach to discretize the lateral skew flux
 \eqref{eq:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995
+\eqref{eq:triad:i13c} from tracer cell $i,k$ to $i+1,k$, introduced in 1995
 into OPA, \eqref{Eq_tra_ldf_iso}, is to calculate a mean vertical
 gradient at the u-point from the average of the four surrounding
+gradient at the $u$-point from the average of the four surrounding
 vertical tracer gradients, and multiply this by a mean slope at the
+u-point, calculated from the averaged surrounding vertical density
+gradients. The total area-integrated skew-flux from tracer cell $i,k$
+to $i+1,k$, noting that the $e_{{3u}_{i+1/2}^k}$ in the area
+$e_{{3u}_{i+1/2}^k}e_{{2u}_{i+1/2}^k}$ at the u-point cancels out with
+the $1/e_{{3u}_{i+1/2}^k}$ associated with the vertical tracer
+$u$-point, calculated from the averaged surrounding vertical density
+gradients. The total area-integrated skew-flux (flux per unit area in
+$ijk$ space) from tracer cell $i,k$
+to $i+1,k$, noting that the $e_{{3}_{i+1/2}^k}$ in the area
+$e{_{3}}_{i+1/2}^k{e_{2}}_{i+1/2}i^k$ at the $u$-point cancels out with
+the $1/{e_{3}}_{i+1/2}^k$ associated with the vertical tracer
 gradient, is then \eqref{Eq_tra_ldf_iso}
 \begin{equation*}
   \left(F_u^{\mathrm{skew}} \right)_{i+\hhalf}^k = \kappa _{i+\hhalf}^k
   e_{{2u}_{i+1/2}^k} \overline{\overline
+  \left(F_u^{13} \right)_{i+\hhalf}^k = \Alts_{i+\hhalf}^k
+  {e_{2}}_{i+1/2}^k \overline{\overline
     r_1} ^{\,i,k}\,\overline{\overline{\delta_k T}}^{\,i,k},
 \end{equation*}
 …
 \begin{equation*}
   \overline{\overline
     r_1} ^{\,i,k} = -\frac{e_{{3u}_{i+1/2}^k}}{e_{{1u}_{i+1/2}^k}}
   \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}}.
+   r_1} ^{\,i,k} = -\frac{{e_{3u}}_{i+1/2}^k}{{e_{1u}}_{i+1/2}^k}
+  \frac{\delta_{i+1/2} [\rho]}{\overline{\overline{\delta_k \rho}}^{\,i,k}},
 \end{equation*}
+and here and in the following we drop the $^{lT}$ superscript from
+$\Alt$ for simplicity.
 Unfortunately the resulting combination $\overline{\overline{\delta_k
     \bullet}}^{\,i,k}$ of a $k$ average and a $k$ difference %of the tracer
 …
 global-average variance. To correct this, we introduced a smoothing of
 the slopes of the iso-neutral surfaces (see \S\ref{LDF}). This
 technique works fine for $T$ and $S$ as they are active tracers
+technique works for $T$ and $S$ in so far as they are active tracers
 ($i.e.$ they enter the computation of density), but it does not work
 for a passive tracer.
 …
 \citep{Griffies_al_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
 % the mean vertical gradient at the u-point,
+% Instead of multiplying the mean slope calculated at the $u$-point by
+% the mean vertical gradient at the $u$-point,
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h] \begin{center}
     \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_triad_fluxes}
     \caption{  \label{Fig_ISO_triad}
+    \caption{ \label{fig:triad:ISO_triad}
       (a) Arrangement of triads $S_i$ and tracer gradients to
            give lateral tracer flux from box $i,k$ to $i+1,k$
+           give lateral tracer flux from box $i,k$ to $i+1,k$
       (b) Triads $S'_i$ and tracer gradients to give vertical tracer flux from
             box $i,k$ to $i,k+1$.}
+    \label{Fig_triad}
+  \end{center} \end{figure}
+ \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 They get the skew flux from the products of the vertical gradients at
 each w-point surrounding the u-point with the corresponding `triad'
 slope calculated from the lateral density gradient across the u-point
 divided by the vertical density gradient at the same w-point as the
 tracer gradient. See Fig.~\ref{Fig_triad}a, where the thick lines
+each $w$-point surrounding the $u$-point with the corresponding `triad'
+slope calculated from the lateral density gradient across the $u$-point
+divided by the vertical density gradient at the same $w$-point as the
+tracer gradient. See Fig.~\ref{fig:triad:ISO_triad}a, where the thick lines
 denote the tracer gradients, and the thin lines the corresponding
 triads, with slopes $s_1, \dotsc s_4$. The total area-integrated
 skew-flux from tracer cell $i,k$ to $i+1,k$
 \begin{multline}
   \label{eq:i13}
   \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \kappa _{i+1}^k a_1 s_1
+  \label{eq:triad:i13}
+  \left( F_u^{13}  \right)_{i+\frac{1}{2}}^k = \Alts_{i+1}^k a_1 s_1
   \delta _{k+\frac{1}{2}} \left[ T^{i+1}
   \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  + \kappa _i^k a_2 s_2 \delta
+  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  + \Alts _i^k a_2 s_2 \delta
   _{k+\frac{1}{2}} \left[ T^i
   \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}} \\
    +\kappa _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1}
   \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  +\kappa _i^k a_4 s_4 \delta
+   +\Alts _{i+1}^k a_3 s_3 \delta _{k-\frac{1}{2}} \left[ T^{i+1}
+  \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}}  +\Alts _i^k a_4 s_4 \delta
   _{k-\frac{1}{2}} \left[ T^i \right]/e_{{3w}_{i+1}}^{k+\frac{1}{2}},
 \end{multline}
 where the contributions of the triad fluxes are weighted by areas
 $a_1, \dotsc a_4$, and $\kappa$ is now defined at the tracer points
 rather than the u-points. This discretization gives a much closer
+$a_1, \dotsc a_4$, and $\Alts$ is now defined at the tracer points
+rather than the $u$-points. This discretization gives a much closer
 stencil, and disallows the two-point computational modes.
  The vertical skew flux \eqref{eq:i31c} from tracer cell $i,k$ to $i,k+1$ at the
 w-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{Fig_triad}b)
+ The vertical skew flux \eqref{eq:triad:i31c} from tracer cell $i,k$ to $i,k+1$ at the
+$w$-point $i,k+\hhalf$ is constructed similarly (Fig.~\ref{fig:triad:ISO_triad}b)
 by multiplying lateral tracer gradients from each of the four
 surrounding u-points by the appropriate triad slope:
+surrounding $u$-points by the appropriate triad slope:
 \begin{multline}
   \label{eq:i31}
   \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \kappa_i^{k+1} a_{1}'
+  \label{eq:triad:i31}
+  \left( F_w^{31} \right) _i ^{k+\frac{1}{2}} =  \Alts_i^{k+1} a_{1}'
   s_{1}' \delta _{i-\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i-\frac{1}{2}}^{k+1}
    +\kappa_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\
   + \kappa_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k
   +\kappa_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k.
+   +\Alts_i^{k+1} a_{2}' s_{2}' \delta _{i+\frac{1}{2}} \left[ T^{k+1} \right]/{e_{3u}}_{i+\frac{1}{2}}^{k+1}\\
+  + \Alts_i^k a_{3}' s_{3}' \delta _{i-\frac{1}{2}} \left[ T^k\right]/{e_{3u}}_{i-\frac{1}{2}}^k
+  +\Alts_i^k a_{4}' s_{4}' \delta _{i+\frac{1}{2}} \left[ T^k \right]/{e_{3u}}_{i+\frac{1}{2}}^k.
 \end{multline}
 We notate the triad slopes in terms of the `anchor point' $i,k$
 (appearing in both the vertical and lateral gradient), and the u- and
 w-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the
 triad as follows (see also Fig.~\ref{Fig_triad}):
 \begin{equation}
   \label{Gf_slopes}
   _i^k \mathbb{R}_{i_p}^{k_p}
   =\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}
+We notate the triad slopes $s_i$ and $s'_i$ in terms of the `anchor point' $i,k$
+(appearing in both the vertical and lateral gradient), and the $u$- and
+$w$-points $(i+i_p,k)$, $(i,k+k_p)$ at the centres of the `arms' of the
+triad as follows (see also Fig.~\ref{fig:triad:ISO_triad}):
+\begin{equation}
+  \label{eq:triad:R}
+  _i^k \mathbb{R}_{i_p}^{k_p}
+  =-\frac{ {e_{3w}}_{\,i}^{\,k+k_p}} { {e_{1u}}_{\,i+i_p}^{\,k}}
+  \
   \frac
+  \frac
   {\left(\alpha / \beta \right)_i^k  \ \delta_{i + i_p}[T^k] - \delta_{i + i_p}[S^k] }
   {\left(\alpha / \beta \right)_i^k  \ \delta_{k+k_p}[T^i ] - \delta_{k+k_p}[S^i ] }.
 …
 In calculating the slopes of the local neutral
 surfaces, the expansion coefficients $\alpha$ and $\beta$ are
 evaluated at the anchor points of the triad \footnote{Note that in \eqref{Gf_slopes} we use the ratio $\alpha / \beta$
+evaluated at the anchor points of the triad \footnote{Note that in \eqref{eq:triad:R} we use the ratio $\alpha / \beta$
 instead of multiplying the temperature derivative by $\alpha$ and the
 salinity derivative by $\beta$. This is more efficient as the ratio
 $\alpha / \beta$ can to be evaluated directly}, while the metrics are
 calculated at the u- and w-points on the arms.
+calculated at the $u$- and $w$-points on the arms.
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h] \begin{center}
     \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_qcells}
+    \caption{   \label{Fig_ISO_triad_notation}
+    Triad notation for quarter cells.T-cells are inside
+      boxes, while the  $i+\half,k$ u-cell is shaded in green and the
+      $i,k+\half$ w-cell is shaded in pink.}
+    \label{qcells}
+    \caption{   \label{fig:triad:qcells}
+    Triad notation for quarter cells.$T$-cells are inside
+      boxes, while the  $i+\half,k$ $u$-cell is shaded in green and the
+      $i,k+\half$ $w$-cell is shaded in pink.}
   \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{qcells}) with the quarter
 cell that is the intersection of the $i,k$ T-cell, the $i+i_p,k$
 u-cell and the $i,k+k_p$ w-cell. Expressing the slopes $s_i$ and
 $s'_i$ in \eqref{eq:i13} and \eqref{eq:i31} in this notation, we have
+Each triad $\{_i^k\:_{i_p}^{k_p}\}$ is associated (Fig.~\ref{fig:triad:qcells}) with the quarter
+cell that is the intersection of the $i,k$ $T$-cell, the $i+i_p,k$
+$u$-cell and the $i,k+k_p$ $w$-cell. Expressing the slopes $s_i$ and
+$s'_i$ in \eqref{eq:triad:i13} and \eqref{eq:triad:i31} in this notation, we have
 e.g.\ $s_1=s'_1={\:}_i^k \mathbb{R}_{1/2}^{1/2}$. Each triad slope $_i^k
 \mathbb{R}_{i_p}^{k_p}$ is used once (as an $s$) to calculate the
 lateral flux along its u-arm, at $(i+i_p,k)$, and then again as an
 $s'$ to calculate the vertical flux along its w-arm at
+lateral flux along its $u$-arm, at $(i+i_p,k)$, and then again as an
+$s'$ to calculate the vertical flux along its $w$-arm at
 $(i,k+k_p)$. Each vertical area $a_i$ used to calculate the lateral
 flux and horizontal area $a'_i$ used to calculate the vertical flux
 can also be identified as the area across the u- and w-arms of a
 unique triad, and we can notate these areas, similarly to the triad
+can also be identified as the area across the $u$- and $w$-arms of a
+unique triad, and we notate these areas, similarly to the triad
 slopes, as $_i^k{\mathbb{A}_u}_{i_p}^{k_p}$,
 $_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:i13}
 $a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:i31}
+$_i^k{\mathbb{A}_w}_{i_p}^{k_p}$, where e.g. in \eqref{eq:triad:i13}
+$a_{1}={\:}_i^k{\mathbb{A}_u}_{1/2}^{1/2}$, and in \eqref{eq:triad:i31}
 $a'_{1}={\:}_i^k{\mathbb{A}_w}_{1/2}^{1/2}$.
 \subsection{The full triad fluxes}
 A key property of isoneutral diffusion is that it should not affect
+A key property of iso-neutral diffusion is that it should not affect
 the (locally referenced) density. In particular there should be no
 lateral or vertical density flux. The lateral density flux disappears so long as the
 …
 form
 \begin{equation}
   \label{eq:i11}
+  \label{eq:triad:i11}
   \left( F_u^{11} \right) _{i+\frac{1}{2}} ^{k} =
   - \left( \kappa_i^{k+1} a_{1} + \kappa_i^{k+1} a_{2} + \kappa_i^k
     a_{3} + \kappa_i^k a_{4} \right)
+  - \left( \Alts_i^{k+1} a_{1} + \Alts_i^{k+1} a_{2} + \Alts_i^k
+    a_{3} + \Alts_i^k a_{4} \right)
   \frac{\delta _{i+1/2} \left[ T^k\right]}{{e_{1u}}_{\,i+1/2}^{\,k}},
 \end{equation}
 where the areas $a_i$ are as in \eqref{eq:i13}. In this case,
 separating the total lateral flux, the sum of \eqref{eq:i13} and
 \eqref{eq:i11}, into triad components, a lateral tracer
+where the areas $a_i$ are as in \eqref{eq:triad:i13}. In this case,
+separating the total lateral flux, the sum of \eqref{eq:triad:i13} and
+\eqref{eq:triad:i11}, into triad components, a lateral tracer
 flux
 \begin{equation}
   \label{eq:latflux-triad}
   _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = -A_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}
+  \label{eq:triad:latflux-triad}
+  _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) = - \Alts_i^k{ \:}_i^k{\mathbb{A}_u}_{i_p}^{k_p}
   \left(
     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
 …
 density flux associated with each triad separately disappears.
 \begin{equation}
   \label{eq:latflux-rho}
+  \label{eq:triad:latflux-rho}
   {\mathbb{F}_u}_{i_p}^{k_p} (\rho)=-\alpha _i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{F}_u}_{i_p}^{k_p} (S)=0
 \end{equation}
 …
 to $i+1,k$ must also vanish since it is a sum of four such triad fluxes.
 The squared slope $r_1^2$ in the expression \eqref{eq:i33c} for the
+The squared slope $r_1^2$ in the expression \eqref{eq:triad:i33c} for the
 $_{33}$ component is also expressed in terms of area-weighted
 squared triad slopes, so the area-integrated vertical flux from tracer
 cell $i,k$ to $i,k+1$ resulting from the $r_1^2$ term is
 \begin{equation}
   \label{eq:i33}
+  \label{eq:triad:i33}
   \left( F_w^{33} \right) _i^{k+\frac{1}{2}} =
     - \left( \kappa_i^{k+1} a_{1}' s_{1}'^2
     + \kappa_i^{k+1} a_{2}' s_{2}'^2
     + \kappa_i^k a_{3}' s_{3}'^2
     + \kappa_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
+    - \left( \Alts_i^{k+1} a_{1}' s_{1}'^2
+    + \Alts_i^{k+1} a_{2}' s_{2}'^2
+    + \Alts_i^k a_{3}' s_{3}'^2
+    + \Alts_i^k a_{4}' s_{4}'^2 \right)\delta_{k+\frac{1}{2}} \left[ T^{i+1} \right],
 \end{equation}
 where the areas $a'$ and slopes $s'$ are the same as in
 \eqref{eq:i31}.
 Then, separating the total vertical flux, the sum of \eqref{eq:i31} and
 \eqref{eq:i33}, into triad components,  a vertical flux
+\eqref{eq:triad:i31}.
+Then, separating the total vertical flux, the sum of \eqref{eq:triad:i31} and
+\eqref{eq:triad:i33}, into triad components,  a vertical flux
 \begin{align}
   \label{eq:vertflux-triad}
+  \label{eq:triad:vertflux-triad}
   _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
   &= A_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}
+  &= \Alts_i^k{\: }_i^k{\mathbb{A}_w}_{i_p}^{k_p}
   \left(
     {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
 …
   \right) \\
   &= - \left(\left.{ }_i^k{\mathbb{A}_w}_{i_p}^{k_p}\right/{ }_i^k{\mathbb{A}_u}_{i_p}^{k_p}\right)
    {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:vertflux-triad2}
+   {_i^k\mathbb{R}_{i_p}^{k_p}}{\: }_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \label{eq:triad:vertflux-triad2}
 \end{align}
 may be associated with each triad. Each vertical density flux $_i^k {\mathbb{F}_w}_{i_p}^{k_p} (\rho)$
 …
 fluxes.
 We can explicitly identify (Fig.~\ref{qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of
 the u-fluxes and w-fluxes in
 \eqref{eq:i31}, \eqref{eq:i13}, \eqref{eq:i11} \eqref{eq:i33} and
 Fig.~\ref{Fig_triad} to  write out the iso-neutral fluxes at $u$- and
+We can explicitly identify (Fig.~\ref{fig:triad:qcells}) the triads associated with the $s_i$, $a_i$, and $s'_i$, $a'_i$ used in the definition of
+the $u$-fluxes and $w$-fluxes in
+\eqref{eq:triad:i31}, \eqref{eq:triad:i13}, \eqref{eq:triad:i11} \eqref{eq:triad:i33} and
+Fig.~\ref{fig:triad:ISO_triad} to  write out the iso-neutral fluxes at $u$- and
 $w$-points as sums of the triad fluxes that cross the $u$- and $w$-faces:
 %(Fig.~\ref{Fig_ISO_triad}):
 \begin{flalign} \label{Eq_iso_flux} \textbf{F}_{iso}(T) &\equiv
+\begin{flalign} \label{Eq_iso_flux} \vect{F}_\mathrm{iso}(T) &\equiv
   \sum_{\substack{i_p,\,k_p}}
   \begin{pmatrix}
 …
   \end{pmatrix}.
 \end{flalign}
 \subsection{Ensuring the scheme cannot increase tracer variance}
 \label{sec:variance}
 We now require that this operator cannot increase the
+\subsection{Ensuring the scheme does not increase tracer variance}
+\label{sec:triad:variance}
+We now require that this operator should not increase the
 globally-integrated tracer variance.
 %This changes according to
 % \begin{align*}
 % &\int_D  D_l^T \; T \;dv \equiv  \sum_{i,k} \left\{ T \ D_l^T \ b_T \right\}    \\
 % &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
 %     \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
+% &\equiv + \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+%     \delta_{i} \left[{_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \right]
 %       + \delta_{k} \left[ {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}} \right]  \ T \right\}    \\
 % &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
+% &\equiv  - \sum_{i,k} \sum_{\substack{i_p,\,k_p}} \left\{
 %                 {_{i+1/2-i_p}^k {\mathbb{F}_u }_{i_p}^{k_p}} \ \delta_{i+1/2} [T]
 %              + {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p}}  \ \delta_{k+1/2} [T]   \right\}      \\
 % \end{align*}
 Each triad slope $_i^k\mathbb{R}_{i_p}^{k_p}$ drives a lateral flux
 $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the u-point $i+i_p,k$ and
+$_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ across the $u$-point $i+i_p,k$ and
 a vertical flux $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ across the
 w-point $i,k+k_p$.  The lateral flux drives a net rate of change of
 variance at points $i+i_p-\half,k$ and $i+i_p+\half,k$ of
+$w$-point $i,k+k_p$.  The lateral flux drives a net rate of change of
+variance, summed over the two $T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
 \begin{multline}
   {b_T}_{i+i_p-1/2}^k\left(\frac{\partial T}{\partial t}T\right)_{i+i_p-1/2}^k+
 …
   &= -T_{i+i_p-1/2}^k{\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \quad + \quad  T_{i+i_p+1/2}^k
   {\;}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T) \\
   &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{locFdtdx}
+  &={\;} _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k], \label{eq:triad:dvar_iso_i}
  \end{split}
 \end{multline}
+while the vertical flux similarly drives a net rate of change of variance at points $i,k+k_p-\half$ and
+$i,k+k_p+\half$ of
+\begin{equation}
+\label{locFdtdz}
+while the vertical flux similarly drives a net rate of change of
+variance summed over the $T$-points $i,k+k_p-\half$ (above) and
+$i,k+k_p+\half$ (below) of
+\begin{equation}
+\label{eq:triad:dvar_iso_k}
   _i^k{\mathbb{F}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
 \end{equation}
 The total variance tendency driven by the triad is the sum of these
 two. Expanding $_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)$ and
 $_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:latflux-triad} and
 \eqref{eq:vertflux-triad}, it is
+$_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)$ with \eqref{eq:triad:latflux-triad} and
+\eqref{eq:triad:vertflux-triad}, it is
 \begin{multline*}
 -A_i^k\left \{
+-\Alts_i^k\left \{
 { } _i^k{\mathbb{A}_u}_{i_p}^{k_p}
   \left(
 …
 to be related to a triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$ by
 \begin{equation}
   \label{eq:V-A}
+  \label{eq:triad:V-A}
   _i^k\mathbb{V}_{i_p}^{k_p}
   ={\;}_i^k{\mathbb{A}_u}_{i_p}^{k_p}\,{e_{1u}}_{\,i + i_p}^{\,k}
 …
 the variance tendency reduces to the perfect square
 \begin{equation}
   \label{eq:perfect-square}
   -A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
+  \label{eq:triad:perfect-square}
+  -\Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
   \left(
     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
 …
   \right)^2\leq 0.
 \end{equation}
 Thus, the constraint \eqref{eq:V-A} ensures that the fluxes associated
+Thus, the constraint \eqref{eq:triad:V-A} ensures that the fluxes (\ref{eq:triad:latflux-triad}, \ref{eq:triad:vertflux-triad}) associated
 with a given slope triad $_i^k\mathbb{R}_{i_p}^{k_p}$ do not increase
 the net variance. Since the total fluxes are sums of such fluxes from
 …
 increase.
 The expression \eqref{eq:V-A} can be interpreted as a discretization
+The expression \eqref{eq:triad:V-A} can be interpreted as a discretization
 of the global integral
 \begin{equation}
   \label{eq:cts-var}
   \frac{\partial}{\partial t}\int\half T^2\, dV =
   \int\mathbf{F}\cdot\nabla T\, dV,
+  \label{eq:triad:cts-var}
+  \frac{\partial}{\partial t}\int\!\half T^2\, dV =
+  \int\!\mathbf{F}\cdot\nabla T\, dV,
 \end{equation}
 where, within each triad volume $_i^k\mathbb{V}_{i_p}^{k_p}$, the
 …
 \[
 \mathbf{F}=\left(
 _i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_u}_{i_p}^{k_p},
 {\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)/_i^k{\mathbb{A}_w}_{i_p}^{k_p}
+\left.{}_i^k{\mathbb{F}_u}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_u}_{i_p}^{k_p},
+\left.{\:}_i^k{\mathbb{F}_w}_{i_p}^{k_p} (T)\right/{}_i^k{\mathbb{A}_w}_{i_p}^{k_p}
  \right)
 \]
 and the gradient
  \[\nabla T = \left(
 \delta_{i+ i_p}[T^k] / {e_{1u}}_{\,i + i_p}^{\,k},
 \delta_{k+ k_p}[T^i] / {e_{3w}}_{\,i}^{\,k + k_p}
+\left.\delta_{i+ i_p}[T^k] \right/ {e_{1u}}_{\,i + i_p}^{\,k},
+\left.\delta_{k+ k_p}[T^i] \right/ {e_{3w}}_{\,i}^{\,k + k_p}
 \right)
 \]
 …
 volumes $_i^k\mathbb{V}_{i_p}^{k_p}$. \citet{Griffies_al_JPO98} identify
 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
 \eqref{eq:cts-var}. The \NEMO model, however, operates with scale
+cells, defined in terms of the distances between $T$, $u$,$f$ and
+$w$-points. This is the natural discretization of
+\eqref{eq:triad:cts-var}. The \NEMO model, however, operates with scale
 factors instead of grid sizes, and scale factors for the quarter
 cells are not defined. Instead, therefore we simply choose
 \begin{equation}
   \label{eq:V-NEMO}
+  \label{eq:triad:V-NEMO}
   _i^k\mathbb{V}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k,
 \end{equation}
 as a quarter of the volume of the u-cell inside which the triad
+as a quarter of the volume of the $u$-cell inside which the triad
 quarter-cell lies. This has the nice property that when the slopes
 $\mathbb{R}$ vanish, the lateral flux from tracer cell $i,k$ to
 $i+1,k$ reduces to the classical form
 \begin{equation}
   \label{eq:lat-normal}
 -\overline{A}_{\,i+1/2}^k\;
+  \label{eq:triad:lat-normal}
+-\overline\Alts_{\,i+1/2}^k\;
 \frac{{b_u}_{i+1/2}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
 \;\frac{\delta_{i+ 1/2}[T^k] }{{e_{1u}}_{\,i + i_p}^{\,k}}
  = -\overline{A}_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}}{{e_{1u}}_{\,i + 1/2}^{\,k}}.
 \end{equation}
 In fact if the diffusive coefficient is defined at u-points, so that
 we employ $A_{i+i_p}^k$ instead of $A_i^k$ in the definitions of the
 triad fluxes \eqref{eq:latflux-triad} and \eqref{eq:vertflux-triad},
+ = -\overline\Alts_{\,i+1/2}^k\;\frac{{e_{1w}}_{\,i + 1/2}^{\,k}\:{e_{1v}}_{\,i + 1/2}^{\,k}\;\delta_{i+ 1/2}[T^k]}{{e_{1u}}_{\,i + 1/2}^{\,k}}.
+\end{equation}
+In fact if the diffusive coefficient is defined at $u$-points, so that
+we employ $\Alts_{i+i_p}^k$ instead of  $\Alts_i^k$ in the definitions of the
+triad fluxes \eqref{eq:triad:latflux-triad} and \eqref{eq:triad:vertflux-triad},
 we can replace $\overline{A}_{\,i+1/2}^k$ by $A_{i+1/2}^k$ in the above.
 \subsection{Summary of the scheme}
+The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
+The iso-neutral fluxes at $u$- and
+$w$-points are the sums of the triad fluxes that cross the $u$- and
+$w$-faces \eqref{Eq_iso_flux}:
+\begin{subequations}\label{eq:triad:alltriadflux}
+  \begin{flalign}\label{eq:triad:vect_isoflux}
+    \vect{F}_\mathrm{iso}(T) &\equiv
+    \sum_{\substack{i_p,\,k_p}}
+    \begin{pmatrix}
+      {_{i+1/2-i_p}^k {\mathbb{F}_u}_{i_p}^{k_p} } (T)      \\
+      \\
+      {_i^{k+1/2-k_p} {\mathbb{F}_w}_{i_p}^{k_p} } (T)
+    \end{pmatrix},
+  \end{flalign}
+  where \eqref{eq:triad:latflux-triad}:
+  \begin{align}
+    \label{eq:triad:triadfluxu}
+    _i^k {\mathbb{F}_u}_{i_p}^{k_p} (T) &= - \Alts_i^k{
+      \:}\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{1u}}_{\,i + i_p}^{\,k}}
+    \left(
+      \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
+      -\ {_i^k\mathbb{R}_{i_p}^{k_p}} \
+      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
+    \right),\\
+    \intertext{and}
+    _i^k {\mathbb{F}_w}_{i_p}^{k_p} (T)
+    &= \Alts_i^k{\: }\frac{{{}_i^k\mathbb{V}}_{i_p}^{k_p}}{{e_{3w}}_{\,i}^{\,k+k_p}}
+    \left(
+      {_i^k\mathbb{R}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
+      -\ \left({_i^k\mathbb{R}_{i_p}^{k_p}}\right)^2 \
+      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} }
+    \right),\label{eq:triad:triadfluxw}
+  \end{align}
+  with \eqref{eq:triad:V-NEMO}
+  \begin{equation}
+    \label{eq:triad:V-NEMO2}
+    _i^k{\mathbb{V}}_{i_p}^{k_p}=\quarter {b_u}_{i+i_p}^k.
+  \end{equation}
+\end{subequations}
+ The divergence of the expression \eqref{Eq_iso_flux} for the fluxes gives the iso-neutral diffusion tendency at
 each tracer point:
 \begin{equation} \label{Gf_operator} D_l^T = \frac{1}{b_T}
+\begin{equation} \label{eq:triad:iso_operator} D_l^T = \frac{1}{b_T}
   \sum_{\substack{i_p,\,k_p}} \left\{ \delta_{i} \left[{_{i+1/2-i_p}^k
         {\mathbb{F}_u }_{i_p}^{k_p}} \right] + \delta_{k} \left[
 …
 \begin{description}
 \item[$\bullet$ horizontal diffusion] The discretization of the
   diffusion operator recovers \eqref{eq:lat-normal} the traditional five-point Laplacian in
+  diffusion operator recovers \eqref{eq:triad:lat-normal} the traditional five-point Laplacian in
   the limit of flat iso-neutral direction :
   \begin{equation} \label{Gf_property1a} D_l^T = \frac{1}{b_T} \
+  \begin{equation} \label{eq:triad:iso_property0} D_l^T = \frac{1}{b_T} \
     \delta_{i} \left[ \frac{e_{2u}\,e_{3u}}{e_{1u}} \;
       \overline{A}^{\,i} \; \delta_{i+1/2}[T] \right] \qquad
+      \overline\Alts^{\,i} \; \delta_{i+1/2}[T] \right] \qquad
     \text{when} \quad { _i^k \mathbb{R}_{i_p}^{k_p} }=0
   \end{equation}
 …
 \item[$\bullet$ implicit treatment in the vertical] Only tracer values
   associated with a single water column appear in the expression
   \eqref{eq:i33} for the $_{33}$ fluxes, vertical fluxes driven by
+  \eqref{eq:triad:i33} for the $_{33}$ fluxes, vertical fluxes driven by
   vertical gradients. This is of paramount importance since it means
+  that an implicit in time algorithm can be used to solve the vertical
+  diffusion equation. This is a necessity since the vertical eddy
+  that a time-implicit algorithm can be used to solve the vertical
+  diffusion equation. This is necessary
+ since the vertical eddy
   diffusivity associated with this term,
   \begin{equation}
     \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{
       {\:}_i^k\mathbb{V}_{i_p}^{k_p} \:A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
     \right\}  =
     \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{
       {b_u}_{i+i_p}^k\:A_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
+    \frac{1}{b_w}\sum_{\substack{i_p, \,k_p}} \left\{
+      {\:}_i^k\mathbb{V}_{i_p}^{k_p} \: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
+    \right\}  =
+    \frac{1}{4b_w}\sum_{\substack{i_p, \,k_p}} \left\{
+      {b_u}_{i+i_p}^k\: \Alts_i^k \: \left(_i^k \mathbb{R}_{i_p}^{k_p}\right)^2
     \right\},
  \end{equation}
 …
 \item[$\bullet$ pure iso-neutral operator] The iso-neutral flux of
   locally referenced potential density is zero. See
   \eqref{eq:latflux-rho} and \eqref{eq:vertflux-triad2}.
+  \eqref{eq:triad:latflux-rho} and \eqref{eq:triad:vertflux-triad2}.
 \item[$\bullet$ conservation of tracer] The iso-neutral diffusion
   conserves tracer content, $i.e.$
   \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ D_l^T \
+  \begin{equation} \label{eq:triad:iso_property1} \sum_{i,j,k} \left\{ D_l^T \
       b_T \right\} = 0
   \end{equation}
 …
 \item[$\bullet$ no increase of tracer variance] The iso-neutral diffusion
   does not increase the tracer variance, $i.e.$
   \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ T \ D_l^T
+  \begin{equation} \label{eq:triad:iso_property2} \sum_{i,j,k} \left\{ T \ D_l^T
       \ b_T \right\} \leq 0
   \end{equation}
   The property is demonstrated in
   \ref{sec:variance} above. It is a key property for a diffusion
   term. It means that it is also a dissipation term, $i.e.$ it is a
   diffusion of the square of the quantity on which it is applied.  It
+  \S\ref{sec:triad:variance} above. It is a key property for a diffusion
+  term. It means that it is also a dissipation term, $i.e.$ it
+  dissipates the square of the quantity on which it is applied.  It
   therefore ensures that, when the diffusivity coefficient is large
   enough, the field on which it is applied become free of grid-point
+  enough, the field on which it is applied becomes free of grid-point
   noise.
 \item[$\bullet$ self-adjoint operator] The iso-neutral diffusion
   operator is self-adjoint, $i.e.$
   \begin{equation} \label{Gf_property1} \sum_{i,j,k} \left\{ S \ D_l^T
+  \begin{equation} \label{eq:triad:iso_property3} \sum_{i,j,k} \left\{ S \ D_l^T
       \ b_T \right\} = \sum_{i,j,k} \left\{ D_l^S \ T \ b_T \right\}
   \end{equation}
 …
   routine. This property can be demonstrated similarly to the proof of
   the `no increase of tracer variance' property. The contribution by a
   single triad towards the left hand side of \eqref{Gf_property1}, can
   be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{locFdtdx}
   and \eqref{locFdtdx}. This results in a term similar to
   \eqref{eq:perfect-square},
 \begin{equation}
   \label{eq:TScovar}
   -A_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
+  single triad towards the left hand side of \eqref{eq:triad:iso_property3}, can
+  be found by replacing $\delta[T]$ by $\delta[S]$ in \eqref{eq:triad:dvar_iso_i}
+  and \eqref{eq:triad:dvar_iso_k}. This results in a term similar to
+  \eqref{eq:triad:perfect-square},
+\begin{equation}
+  \label{eq:triad:TScovar}
+  - \Alts_i^k{\:} _i^k\mathbb{V}_{i_p}^{k_p}
   \left(
     \frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }
 …
 This is symmetrical in $T $ and $S$, so exactly the same term arises
 from the discretization of this triad's contribution towards the
 RHS of \eqref{Gf_property1}.
+RHS of \eqref{eq:triad:iso_property3}.
 \end{description}
+$\ $\newline %force an empty line
+\subsection{Treatment of the triads at the boundaries}\label{sec:triad:iso_bdry}
+The triad slope can only be defined where both the grid boxes centred at
+the end of the arms exist. Triads that would poke up
+through the upper ocean surface into the atmosphere, or down into the
+ocean floor, must be masked out. See Fig.~\ref{fig:triad:bdry_triads}. Surface layer triads
+$\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and
+$\triad{i+1}{1}{R}{-1/2}{-1/2}$ (blue) that require density to be
+specified above the ocean surface are masked (Fig.~\ref{fig:triad:bdry_triads}a): this ensures that lateral
+tracer gradients produce no flux through the ocean surface. However,
+to prevent surface noise, it is customary to retain the $_{11}$ contributions towards
+the lateral triad fluxes $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and
+$\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$; this drives diapycnal tracer
+fluxes. Similar comments apply to triads that would intersect the
+ocean floor (Fig.~\ref{fig:triad:bdry_triads}b). Note that both near bottom
+triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
+$\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
+or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point is
+masked. The associated lateral fluxes (grey-black dashed line) are
+masked if \np{ln\_botmix\_grif}=false, but left unmasked,
+giving bottom mixing, if \np{ln\_botmix\_grif}=true.
+The default option \np{ln\_botmix\_grif}=false is suitable when the
+bbl mixing option is enabled (\key{trabbl}, with \np{nn\_bbl\_ldf}=1),
+or  for simple idealized  problems. For setups with topography without
+bbl mixing, \np{ln\_botmix\_grif}=true may be necessary.
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[h] \begin{center}
+    \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_bdry_triads}
+    \caption{  \label{fig:triad:bdry_triads}
+      (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer
+      points (black dots), and $i+1/2,1$ $u$-point (blue square). Triad
+      slopes $\triad{i}{1}{R}{1/2}{-1/2}$ (magenta) and $\triad{i+1}{1}{R}{-1/2}{-1/2}$
+      (blue) poking through the ocean surface are masked (faded in
+      figure). However, the lateral $_{11}$ contributions towards
+      $\triad[u]{i}{1}{F}{1/2}{-1/2}$ and $\triad[u]{i+1}{1}{F}{-1/2}{-1/2}$
+      (yellow line) are still applied, giving diapycnal diffusive
+      fluxes.\\
+      (b) Both near bottom triad slopes $\triad{i}{k}{R}{1/2}{1/2}$ and
+      $\triad{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the $i,k+1$
+      or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$ $u$-point
+      is masked. The associated lateral fluxes (grey-black dashed
+      line) are masked if \np{botmix\_grif}=.false., but left
+      unmasked, giving bottom mixing, if \np{botmix\_grif}=.true.}
+ \end{center} \end{figure}
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\subsection{ Limiting of the slopes within the interior}\label{sec:triad:limit}
+As discussed in \S\ref{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}. 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.
+%, ramping linearly down above 70~m depth to zero at the surface
+It is of course relevant to the iso-neutral slopes $\tilde{r}_i=r_i+\sigma_i$ relative to
+geopotentials (here the $\sigma_i$ are the slopes of the coordinate surfaces relative to
+geopotentials) \eqref{Eq_PE_slopes_eiv} rather than the slope $r_i$ relative to coordinate
+surfaces, so we require
+\begin{equation*}
+  |\tilde{r}_i|\leq \tilde{r}_\mathrm{max}=0.01.
+\end{equation*}
+and then recalculate the slopes $r_i$ relative to coordinates.
+Each individual triad slope
+ \begin{equation}
+   \label{eq:triad:Rtilde}
+_i^k\tilde{\mathbb{R}}_{i_p}^{k_p} = {}_i^k\mathbb{R}_{i_p}^{k_p}  + \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
+ \end{equation}
+is limited like this and then the corresponding
+$_i^k\mathbb{R}_{i_p}^{k_p} $ are recalculated and combined to form the fluxes.
+Note that where the slopes have been limited, there is now a non-zero
+iso-neutral density flux that drives dianeutral mixing.  In particular this iso-neutral density flux
+is always downwards, and so acts to reduce gravitational potential energy.
+\subsection{Tapering within the surface mixed layer}\label{sec:triad:taper}
+Additional tapering of the iso-neutral fluxes is necessary within the
+surface mixed layer. When the Griffies triads are used, we offer two
+options for this.
+\subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper}
+This is the option activated by the default choice
+\np{ln\_triad\_iso}=false. Slopes $\tilde{r}_i$ relative to
+geopotentials are tapered linearly from their value immediately below the mixed layer to zero at the
+surface, as described in option (c) of Fig.~\ref{Fig_eiv_slp}, to values
+\begin{subequations}
+  \begin{equation}
+   \label{eq:triad:rmtilde}
+     \rMLt =
+  -\frac{z}{h}\left.\tilde{r}_i\right|_{z=-h}\quad \text{ for  } z>-h,
+  \end{equation}
+and then the $r_i$ relative to vertical coordinate surfaces are appropriately
+adjusted to
+  \begin{equation}
+   \label{eq:triad:rm}
+ \rML =\rMLt -\sigma_i \quad \text{ for  } z>-h.
+  \end{equation}
+\end{subequations}
+Thus the diffusion operator within the mixed layer is given by:
+\begin{equation} \label{eq:triad:iso_tensor_ML}
+D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
+\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c}
+\hfill & 0 \hfill & {-\rML[1]}\hfill \\
+\hfill & 1 \hfill & {-\rML[2]} \hfill \\
+ {-\rML[1]}\hfill &   {-\rML[2]} \hfill & {\rML[1]^2+\rML[2]^2} \hfill
+\end{array} }} \right)
+\end{equation}
+This slope tapering gives a natural connection between tracer in the
+mixed-layer and in isopycnal layers immediately below, in the
+thermocline. It is consistent with the way the $\tilde{r}_i$ are
+tapered within the mixed layer (see \S\ref{sec:triad:taperskew} below)
+so as to ensure a uniform GM eddy-induced velocity throughout the
+mixed layer. However, it gives a downwards density flux and so acts so
+as to reduce potential energy in the same way as does the slope
+limiting discussed above in \S\ref{sec:triad:limit}.
+As in \S\ref{sec:triad:limit} above, the tapering
+\eqref{eq:triad:rmtilde} is applied separately to each triad
+$_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}$, and the
+$_i^k\mathbb{R}_{i_p}^{k_p}$ adjusted. For clarity, we assume
+$z$-coordinates in the following; the conversion from
+$\mathbb{R}$ to $\tilde{\mathbb{R}}$ and back to $\mathbb{R}$ follows exactly as described
+above by \eqref{eq:triad:Rtilde}.
+\begin{enumerate}
+\item Mixed-layer depth is defined so as to avoid including regions of weak
+vertical stratification in the slope definition.
+ At each $i,j$ (simplified to $i$ in
+Fig.~\ref{fig:triad:MLB_triad}), we define the mixed-layer by setting
+the vertical index of the tracer point immediately below the mixed
+layer, $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point)
+such that the potential density
+${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$, where $i,k_{10}$ is
+the tracer gridbox within which the depth reaches 10~m. See the left
+side of Fig.~\ref{fig:triad:MLB_triad}. We use the $k_{10}$-gridbox
+instead of the surface gridbox to avoid problems e.g.\ with thin
+daytime mixed-layers. Currently we use the same
+$\Delta\rho_c=0.01\;\mathrm{kg\:m^{-3}}$ for ML triad tapering as is
+used to output the diagnosed mixed-layer depth
+$h_\mathrm{ML}=|z_{W}|_{k_\mathrm{ML}+1/2}$, the depth of the $w$-point
+above the $i,k_\mathrm{ML}$ tracer point.
+\item We define `basal' triad slopes
+${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ as the slopes
+of those triads whose vertical `arms' go down from the
+$i,k_\mathrm{ML}$ tracer point to the $i,k_\mathrm{ML}-1$ tracer point
+below. This is to ensure that the vertical density gradients
+associated with these basal triad slopes
+${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ are
+representative of the thermocline. The four basal triads defined in the bottom part
+of Fig.~\ref{fig:triad:MLB_triad} are then
+\begin{align}
+  {\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p} &=
+ {\:}^{k_\mathrm{ML}-k_p-1/2}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}, \label{eq:triad:Rbase}
+\\
+\intertext{with e.g.\ the green triad}
+{\:}_i{\mathbb{R}_\mathrm{base}}_{1/2}^{-1/2}&=
+{\:}^{k_\mathrm{ML}}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2}. \notag
+\end{align}
+The vertical flux associated with each of these triads passes through the $w$-point
+$i,k_\mathrm{ML}-1/2$ lying \emph{below} the $i,k_\mathrm{ML}$ tracer point,
+so it is this depth
+\begin{equation}
+  \label{eq:triad:zbase}
+  {z_\mathrm{base}}_{\,i}={z_{w}}_{k_\mathrm{ML}-1/2}
+\end{equation}
+(one gridbox deeper than the
+diagnosed ML depth $z_\mathrm{ML})$ that sets the $h$ used to taper
+the slopes in \eqref{eq:triad:rmtilde}.
+\item Finally, we calculate the adjusted triads
+${\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p}$ within the mixed
+layer, by multiplying the appropriate
+${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$ by the ratio of
+the depth of the $w$-point ${z_w}_{k+k_p}$ to ${z_\mathrm{base}}_{\,i}$. For
+instance the green triad centred on $i,k$
+\begin{align}
+  {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,1/2}^{-1/2} &=
+\frac{{z_w}_{k-1/2}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,1/2}^{-1/2}
+\notag \\
+\intertext{and more generally}
+ {\:}_i^k{\mathbb{R}_\mathrm{ML}}_{\,i_p}^{k_p} &=
+\frac{{z_w}_{k+k_p}}{{z_\mathrm{base}}_{\,i}}{\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}.\label{eq:triad:RML}
+\end{align}
+\end{enumerate}
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[h]
+  \fcapside {\caption{\label{fig:triad:MLB_triad} Definition of
+      mixed-layer depth and calculation of linearly tapered
+      triads. The figure shows a water column at a given $i,j$
+      (simplified to $i$), with the ocean surface at the top. Tracer points are denoted by
+      bullets, and black lines the edges of the tracer cells; $k$
+      increases upwards. \\
+      \hspace{5 em}We define the mixed-layer by setting the vertical index
+      of the tracer point immediately below the mixed layer,
+      $k_\mathrm{ML}$, as the maximum $k$ (shallowest tracer point)
+      such that ${\rho_0}_{i,k}>{\rho_0}_{i,k_{10}}+\Delta\rho_c$,
+      where $i,k_{10}$ is the tracer gridbox within which the depth
+      reaches 10~m. We calculate the triad slopes within the mixed
+      layer by linearly tapering them from zero (at the surface) to
+      the `basal' slopes, the slopes of the four triads passing through the
+      $w$-point $i,k_\mathrm{ML}-1/2$ (blue square),
+      ${\:}_i{\mathbb{R}_\mathrm{base}}_{\,i_p}^{k_p}$. Triads with
+    different $i_p,k_p$, denoted by different colours, (e.g. the green
+    triad $i_p=1/2,k_p=-1/2$) are tapered to the appropriate basal triad.}}
+  {\includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_triad_MLB}}
+\end{figure}
+% >>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\subsubsection{Additional truncation of skew iso-neutral flux
+  components}
+\label{sec:triad:Gerdes-taper}
+The alternative option is activated by setting \np{ln\_triad\_iso} =
+  true. This retains the same tapered slope $\rML$  described above for the
+calculation of the $_{33}$ term of the iso-neutral diffusion tensor (the
+vertical tracer flux driven by vertical tracer gradients), but
+replaces the $\rML$ in the skew term by
+\begin{equation}
+  \label{eq:triad:rm*}
+  \rML^*=\left.\rMLt^2\right/\tilde{r}_i-\sigma_i,
+\end{equation}
+giving a ML diffusive operator
+\begin{equation} \label{eq:triad:iso_tensor_ML2}
+D^{lT}=\nabla {\rm {\bf .}}\left( {A^{lT}\;\Re \;\nabla T} \right) \qquad
+\mbox{with}\quad \;\;\Re =\left( {{\begin{array}{*{20}c}
+\hfill & 0 \hfill & {-\rML[1]^*}\hfill \\
+\hfill & 1 \hfill & {-\rML[2]^*} \hfill \\
+ {-\rML[1]^*}\hfill &   {-\rML[2]^*} \hfill & {\rML[1]^2+\rML[2]^2} \hfill \\
+\end{array} }} \right).
+\end{equation}
+This operator
+\footnote{To ensure good behaviour where horizontal density
+  gradients are weak, we in fact follow \citet{Gerdes1991} and set
+$\rML^*=\mathrm{sgn}(\tilde{r}_i)\min(|\rMLt^2/\tilde{r}_i|,|\tilde{r}_i|)-\sigma_i$.}
+then has the property it gives no vertical density flux, and so does
+not change the potential energy.
+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}).
+Again it is applied separately to each triad $_i^k\mathbb{R}_{i_p}^{k_p}$
+In practice, this approach gives weak vertical tracer fluxes through
+the mixed-layer, as well as vanishing density fluxes. While it is
+theoretically advantageous that it does not change the potential
+energy, it may give a discontinuity between the
+fluxes within the mixed-layer (purely horizontal) and just below (along
+iso-neutral surfaces).
+% This may give strange looking results,
+% particularly where the mixed-layer depth varies strongly laterally.
 % ================================================================
 % Skew flux formulation for Eddy Induced Velocity :
 % ================================================================
+\section{ Eddy induced velocity and Skew flux formulation}
+When Gent and McWilliams's [1990] diffusion is used (\key{traldf\_eiv} defined),
+an additional advection term is added. The associated velocity is the so called
+\section{Eddy induced advection formulated as a skew flux}\label{sec:triad:skew-flux}
+\subsection{The continuous skew flux formulation}\label{sec:triad:continuous-skew-flux}
+ When Gent and McWilliams's [1990] diffusion is used,
+an additional advection term is added. The associated velocity is the so called
 eddy induced velocity, the formulation of which depends on the slopes of iso-
 neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used
 here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}
+neutral surfaces. Contrary to the case of iso-neutral mixing, the slopes used
+here are referenced to the geopotential surfaces, $i.e.$ \eqref{Eq_ldfslp_geo}
 is used in $z$-coordinate, and the sum \eqref{Eq_ldfslp_geo}
++ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.
+The eddy induced velocity is given by:
+\begin{equation} \label{Eq_eiv_v}
++ \eqref{Eq_ldfslp_iso} in $z^*$ or $s$-coordinates.
+The eddy induced velocity is given by:
+\begin{subequations} \label{eq:triad:eiv}
+\begin{equation}\label{eq:triad:eiv_v}
 \begin{split}
  u^* & = - \frac{1}{e_{2}e_{3}}\;          \partial_k \left( e_{2} \, A_{e} \; r_i \right)   \\
  v^* & = - \frac{1}{e_{1}e_{3}}\;          \partial_k \left( e_{1} \, A_{e} \; r_j \right)   \\
 w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i  \left( e_{2} \, A_{e} \; r_i \right)
                          + \partial_j  \left( e_{1} \, A_{e} \;r_j \right) \right\} \\
+ u^* & = - \frac{1}{e_{3}}\;          \partial_i\psi_1,  \\
+ v^* & = - \frac{1}{e_{3}}\;          \partial_j\psi_2,    \\
+w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i  \left( e_{2} \, \psi_1\right)
+                         + \partial_j  \left( e_{1} \, \psi_2\right) \right\},
 \end{split}
 \end{equation}
+where $A_{e}$ is the eddy induced velocity coefficient, and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces.
+The traditional way to implement this additional advection is to add it to the Eulerian
+velocity prior to computing the tracer advection. This allows us to take advantage of
+all the advection schemes offered for the tracers (see \S\ref{TRA_adv}) and not just
+a $2^{nd}$ order advection scheme. This is particularly useful for passive tracers
+where \emph{positivity} of the advection scheme is of paramount importance.
+\citet{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
+where the streamfunctions $\psi_i$ are given by
+\begin{equation} \label{eq:triad:eiv_psi}
+\begin{split}
+\psi_1 & = A_{e} \; \tilde{r}_1,   \\
+\psi_2 & = A_{e} \; \tilde{r}_2,
+\end{split}
+\end{equation}
+\end{subequations}
+with $A_{e}$ the eddy induced velocity coefficient, and $\tilde{r}_1$ and $\tilde{r}_2$ the slopes between the iso-neutral and the geopotential surfaces.
+The traditional way to implement this additional advection is to add
+it to the Eulerian velocity prior to computing the tracer
+advection. This is implemented if \key{traldf\_eiv} is set in the
+default implementation, where \np{ln\_traldf\_grif} is set
+false. This allows us to take advantage of all the advection schemes
+offered for the tracers (see \S\ref{TRA_adv}) and not just a $2^{nd}$
+order advection scheme. This is particularly useful for passive
+tracers where \emph{positivity} of the advection scheme is of
+paramount importance.
+However, when \np{ln\_traldf\_grif} is set true, \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, the tracer advective
+fluxes per unit area in $ijk$ space can be
 transformed as follows:
 \begin{flalign*}
 \begin{split}
 \textbf{F}_{eiv}^T =
 \begin{pmatrix}
+\textbf{F}_\mathrm{eiv}^T =
+\begin{pmatrix}
            {e_{2}\,e_{3}\;  u^*}       \\
       {e_{1}\,e_{2}\; w^*}  \\
 \end{pmatrix}   \;   T
 &=
 \begin{pmatrix}
            { - \partial_k \left( e_{2} \, A_{e} \; r_i \right) \; T \;}       \\
       {+ \partial_i  \left( e_{2} \, A_{e} \; r_i \right) \; T \;}    \\
+\begin{pmatrix}
+           { - \partial_k \left( e_{2} \,\psi_1 \right) \; T \;}     \\
+      {+ \partial_i  \left( e_{2} \, \psi_1 \right) \; T \;}    \\
 \end{pmatrix}        \\
 &=
 \begin{pmatrix}
            { - \partial_k \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}  \\
       {+ \partial_i  \left( e_{2} \, A_{e} \; r_i  \; T \right) \;}   \\
 \end{pmatrix}
+ +
 \begin{pmatrix}
            {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}  \\
       { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}  \\
 \end{pmatrix}
+&=
+\begin{pmatrix}
+           { - \partial_k \left( e_{2} \, \psi_1  \; T \right) \;}  \\
+      {+ \partial_i  \left( e_{2} \,\psi_1 \; T \right) \;}  \\
+\end{pmatrix}
+ +
+\begin{pmatrix}
+           {+ e_{2} \, \psi_1  \; \partial_k T}  \\
+      { - e_{2} \, \psi_1  \; \partial_i  T}  \\
+\end{pmatrix}
 \end{split}
 \end{flalign*}
 and since the eddy induces velocity field is no-divergent, we end up with the skew
 form of the eddy induced advective fluxes:
 \begin{equation} \label{Eq_eiv_skew_continuous}
 \textbf{F}_{eiv}^T = \begin{pmatrix}
            {+ e_{2} \, A_{e} \; r_i  \; \partial_k T}   \\
       { - e_{2} \, A_{e} \; r_i  \; \partial_i  T}  \\
+and since the eddy induced velocity field is non-divergent, we end up with the skew
+form of the eddy induced advective fluxes per unit area in $ijk$ space:
+\begin{equation} \label{eq:triad:eiv_skew_ijk}
+\textbf{F}_\mathrm{eiv}^T = \begin{pmatrix}
+           {+ e_{2} \, \psi_1  \; \partial_k T}   \\
+      { - e_{2} \, \psi_1  \; \partial_i  T}  \\
                                  \end{pmatrix}
 \end{equation}
+The tendency associated with eddy induced velocity is then simply the divergence
+of the \eqref{Eq_eiv_skew_continuous} fluxes. It naturally conserves the tracer
+content, as it is expressed in flux form. It also preserve the tracer variance.
+The discrete form of  \eqref{Eq_eiv_skew_continuous} using the slopes \eqref{Gf_slopes} and defining $A_e$ at $T$-point is then given by:
+\begin{flalign} \label{Eq_eiv_skew}   \begin{split}
+\textbf{F}_{eiv}^T   \equiv
+                        \sum_{\substack{i_p,\,k_p}} \begin{pmatrix}
++{e_{2u}}_{i+1/2-i_p}^{k}                                    \ {A_{e}}_{i+1/2-i_p}^{k}
+\ \ { _{i+1/2-i_p}^k \mathbb{R}_{i_p}^{k_p} }    \ \ \delta_{k+k_p}[T_{i+1/2-i_p}]      \\
+    \\
+- {e_{2u}}_i^{k+1/2-k_p}                                      \ {A_{e}}_i^{k+1/2-k_p}
+\ \ { _i^{k+1/2-k_p} \mathbb{R}_{i_p}^{k_p} }    \ \delta_{i+i_p}[T^{k+1/2-k_p}]    \\
+                                                                       \end{pmatrix}
+\end{split}   \end{flalign}
+Note that \eqref{Eq_eiv_skew} is valid in $z$-coordinate with or without partial cells. In $z^*$ or $s$-coordinate, the the slope between the level and the geopotential surfaces must be added to $\mathbb{R}$ for the discret form to be exact.
+Such a choice of discretisation is consistent with the iso-neutral operator as it uses the same definition for the slopes. It also ensures the conservation of the tracer variance (see Appendix \ref{Apdx_eiv_skew}), $i.e.$ it does not include a diffusive component but is a `pure' advection term.
+\newpage      %force an empty line
+% ================================================================
+% Discrete Invariants of the skew flux formulation
+% ================================================================
+\subsection{Discrete Invariants of the skew flux formulation}
+\label{Apdx_eiv_skew}
+Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
+This have to be moved in an Appendix.
+The continuous property to be demonstrated is :
+\begin{align*}
+\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
+\end{align*}
+The discrete form of its left hand side is obtained using \eqref{Eq_eiv_skew}
+\begin{align*}
+ \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
+ \delta_i  &\left[
+{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
+\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
+   \right] \; T_i^k      \\
+- \delta_k &\left[
+{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
+\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
+   \right] \; T_i^k      \         \Biggr\}
+\end{align*}
+apply the adjoint of delta operator, it becomes
+\begin{align*}
+ \sum\limits_{i,k} \sum_{\substack{i_p,\,k_p}}  \Biggl\{   \;\;
+  &\left(
+{e_{2u}}_{i+i_p+1/2}^{k}                                  \;\ \ {A_{e}}_{i+i_p+1/2}^{k}
+\ \ \ { _{i+i_p+1/2}^k \mathbb{R}_{-i_p}^{k_p} }   \quad \delta_{k+k_p}[T_{i+i_p+1/2}]
+   \right) \; \delta_{i+1/2}[T^{k}]      \\
+- &\left(
+{e_{2u}}_i^{k+k_p+1/2}                                     \ \ {A_{e}}_i^{k+k_p+1/2}
+\ \ { _i^{k+k_p+1/2} \mathbb{R}_{i_p}^{-k_p} }   \ \ \delta_{i+i_p}[T^{k+k_p+1/2}]
+     \right) \; \delta_{k+1/2}[T_{i}]       \         \Biggr\}
+\end{align*}
+Expending the summation on $i_p$ and $k_p$, it becomes:
+\begin{align*}
+ \begin{matrix}
+&\sum\limits_{i,k}   \Bigl\{
+  &+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
+  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{-1/2}} &\delta_{k-1/2}[T_{i+1}]    &\delta_{i+1/2}[T^{k}]   &\\
+&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
+  &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{-1/2}}  &\delta_{k-1/2}[T_{i\ \ \ \;}]  &\delta_{i+1/2}[T^{k}] &\\
+&&+{e_{2u}}_{i+1}^{k}                             &{A_{e}}_{i+1    }^{k}
+  &\ {_{i+1}^k \mathbb{R}_{- 1/2}^{+1/2}} &\delta_{k+1/2}[T_{i+1}]     &\delta_{i+1/2}[T^{k}] &\\
+&&+{e_{2u}}_i^{k\ \ \ \:}                            &{A_{e}}_{i}^{k\ \ \ \:}
+    &\ {\ \ \;_i^k \mathbb{R}_{+1/2}^{+1/2}} &\delta_{k+1/2}[T_{i\ \ \ \;}] &\delta_{i+1/2}[T^{k}] &\\
+%
+&&-{e_{2u}}_i^{k+1}                                &{A_{e}}_i^{k+1}
+  &{_i^{k+1} \mathbb{R}_{-1/2}^{- 1/2}}   &\delta_{i-1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
+&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
+  &{\ \ \;_i^k  \mathbb{R}_{-1/2}^{+1/2}}   &\delta_{i-1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}] &\\
+&&-{e_{2u}}_i^{k+1    }                             &{A_{e}}_i^{k+1}
+  &{_i^{k+1} \mathbb{R}_{+1/2}^{- 1/2}}   &\delta_{i+1/2}[T^{k+1}]      &\delta_{k+1/2}[T_{i}] &\\
+&&-{e_{2u}}_i^{k\ \ \ \:}                             &{A_{e}}_i^{k\ \ \ \:}
+  &{\ \ \;_i^k  \mathbb{R}_{+1/2}^{+1/2}}   &\delta_{i+1/2}[T^{k\ \ \ \:}]  &\delta_{k+1/2}[T_{i}]
+&\Bigr\}  \\
+\end{matrix}
+\end{align*}
+The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{+1/2}}$ are the
+same but of opposite signs, they cancel out.
+Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{-1/2}}$.
+The two terms associated with the triad ${_i^k \mathbb{R}_{+1/2}^{-1/2}}$ are the
+same but both of opposite signs and shifted by 1 in $k$ direction. When summing over $k$
+they cancel out with the neighbouring grid points.
+Exactly the same thing occurs for the triad ${_i^k \mathbb{R}_{-1/2}^{+1/2}}$ in the
+$i$ direction. Therefore the sum over the domain is zero, $i.e.$ the variance of the
+tracer is preserved by the discretisation of the skew fluxes.
+%%% Local Variables:
+%%% TeX-master: "../../NEMO_book.tex"
+%%% End:
+The total fluxes per unit physical area are then
+\begin{equation}\label{eq:triad:eiv_skew_physical}
+\begin{split}
+ f^*_1 & = \frac{1}{e_{3}}\; \psi_1 \partial_k T   \\
+ f^*_2 & = \frac{1}{e_{3}}\; \psi_2 \partial_k T   \\
+ f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; \left\{ e_{2} \psi_1 \partial_i T
+   + e_{1} \psi_2 \partial_j T \right\}. \\
+\end{split}
+\end{equation}
+Note that Eq.~ \eqref{eq:triad:eiv_skew_physical} takes the same form whatever the
+vertical coordinate, though of course the slopes
+$\tilde{r}_i$ which define the $\psi_i$ in \eqref{eq:triad:eiv_psi} are relative to geopotentials.
+The tendency associated with eddy induced velocity is then simply the convergence
+of the fluxes (\ref{eq:triad:eiv_skew_ijk}, \ref{eq:triad:eiv_skew_physical}), so
+\begin{equation} \label{eq:triad:skew_eiv_conv}
+\frac{\partial T}{\partial t}= -\frac{1}{e_1 \, e_2 \, e_3 }      \left[
+  \frac{\partial}{\partial i} \left( e_2 \psi_1 \partial_k T\right)
+  + \frac{\partial}{\partial j} \left( e_1  \;
+    \psi_2 \partial_k T\right)
+ -  \frac{\partial}{\partial k} \left( e_{2} \psi_1 \partial_i T
+   + e_{1} \psi_2 \partial_j T \right)  \right]
+\end{equation}
+ It naturally conserves the tracer content, as it is expressed in flux
+ form. Since it has the same divergence as the advective form it also
+ preserves the tracer variance.
+\subsection{The discrete skew flux formulation}
+The skew fluxes in (\ref{eq:triad:eiv_skew_physical}, \ref{eq:triad:eiv_skew_ijk}), like the off-diagonal terms
+(\ref{eq:triad:i13c}, \ref{eq:triad:i31c}) of the small angle diffusion tensor, are best
+expressed in terms of the triad slopes, as in Fig.~\ref{fig:triad:ISO_triad}
+and Eqs~(\ref{eq:triad:i13}, \ref{eq:triad:i31}); but now in terms of the triad slopes
+$\tilde{\mathbb{R}}$ relative to geopotentials instead of the
+$\mathbb{R}$ relative to coordinate surfaces. The discrete form of
+\eqref{eq:triad:eiv_skew_ijk} using the slopes \eqref{eq:triad:R} and
+defining $A_e$ at $T$-points is then given by:
+\begin{subequations}\label{eq:triad:allskewflux}
+  \begin{flalign}\label{eq:triad:vect_skew_flux}
+    \vect{F}_\mathrm{eiv}(T) &\equiv
+    \sum_{\substack{i_p,\,k_p}}
+    \begin{pmatrix}
+      {_{i+1/2-i_p}^k {\mathbb{S}_u}_{i_p}^{k_p} } (T)      \\
+      \\
+      {_i^{k+1/2-k_p} {\mathbb{S}_w}_{i_p}^{k_p} } (T)      \\
+    \end{pmatrix},
+  \end{flalign}
+  where the skew flux in the $i$-direction associated with a given
+  triad is (\ref{eq:triad:latflux-triad}, \ref{eq:triad:triadfluxu}):
+  \begin{align}
+    \label{eq:triad:skewfluxu}
+    _i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) &= + \quarter {A_e}_i^k{
+      \:}\frac{{b_u}_{i+i_p}^k}{{e_{1u}}_{\,i + i_p}^{\,k}}
+     \ {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}} \
+      \frac{ \delta_{k+k_p} [T^i] }{{e_{3w}}_{\,i}^{\,k+k_p} },
+   \\
+    \intertext{and \eqref{eq:triad:triadfluxw} in the $k$-direction, changing the sign
+      to be consistent with \eqref{eq:triad:eiv_skew_ijk}:}
+    _i^k {\mathbb{S}_w}_{i_p}^{k_p} (T)
+    &= -\quarter {A_e}_i^k{\: }\frac{{b_u}_{i+i_p}^k}{{e_{3w}}_{\,i}^{\,k+k_p}}
+     {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}\frac{ \delta_{i+ i_p}[T^k] }{ {e_{1u}}_{\,i + i_p}^{\,k} }.\label{eq:triad:skewfluxw}
+  \end{align}
+\end{subequations}
+Such a discretisation is consistent with the iso-neutral
+operator as it uses the same definition for the slopes.  It also
+ensures the following two key properties.
+\subsubsection{No change in tracer variance}
+The discretization conserves tracer variance, $i.e.$ it does not
+include a diffusive component but is a `pure' advection term. This can
+be seen
+%either from Appendix \ref{Apdx_eiv_skew} or
+by considering the
+fluxes associated with a given triad slope
+$_i^k{\mathbb{R}}_{i_p}^{k_p} (T)$. For, following
+\S\ref{sec:triad:variance} and \eqref{eq:triad:dvar_iso_i}, the
+associated horizontal skew-flux $_i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)$
+drives a net rate of change of variance, summed over the two
+$T$-points $i+i_p-\half,k$ and $i+i_p+\half,k$, of
+\begin{equation}
+\label{eq:triad:dvar_eiv_i}
+  _i^k{\mathbb{S}_u}_{i_p}^{k_p} (T)\,\delta_{i+ i_p}[T^k],
+\end{equation}
+while the associated vertical skew-flux gives a variance change summed over the
+$T$-points $i,k+k_p-\half$ (above) and $i,k+k_p+\half$ (below) of
+\begin{equation}
+\label{eq:triad:dvar_eiv_k}
+  _i^k{\mathbb{S}_w}_{i_p}^{k_p} (T) \,\delta_{k+ k_p}[T^i].
+\end{equation}
+Inspection of the definitions (\ref{eq:triad:skewfluxu}, \ref{eq:triad:skewfluxw})
+shows that these two variance changes (\ref{eq:triad:dvar_eiv_i}, \ref{eq:triad:dvar_eiv_k})
+sum to zero. Hence the two fluxes associated with each triad make no
+net contribution to the variance budget.
+\subsubsection{Reduction in gravitational PE}
+The vertical density flux associated with the vertical skew-flux
+always has the same sign as the vertical density gradient; thus, so
+long as the fluid is stable (the vertical density gradient is
+negative) the vertical density flux is negative (downward) and hence
+reduces the gravitational PE.
+For the change in gravitational PE driven by the $k$-flux is
+\begin{align}
+  \label{eq:triad:vert_densityPE}
+  g {e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho)
+  &=g {e_{3w}}_{\,i}^{\,k+k_p}\left[-\alpha _i^k {\:}_i^k
+    {\mathbb{S}_w}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k
+    {\mathbb{S}_w}_{i_p}^{k_p} (S) \right]. \notag \\
+\intertext{Substituting  ${\:}_i^k {\mathbb{S}_w}_{i_p}^{k_p}$ from
+  \eqref{eq:triad:skewfluxw}, gives}
+% and separating out
+% $\rtriadt{R}=\rtriad{R} + \delta_{i+i_p}[z_T^k]$,
+% gives two terms. The
+% first $\rtriad{R}$ term (the only term for $z$-coordinates) is:
+ &=-\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k {_i^k\tilde{\mathbb{R}}_{i_p}^{k_p}}
+\frac{ -\alpha _i^k\delta_{i+ i_p}[T^k]+ \beta_i^k\delta_{i+ i_p}[S^k]} { {e_{1u}}_{\,i + i_p}^{\,k} } \notag \\
+ &=+\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
+     \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right) {_i^k\mathbb{R}_{i_p}^{k_p}}
+\frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
+\end{align}
+using the definition of the triad slope $\rtriad{R}$,
+\eqref{eq:triad:R} to express $-\alpha _i^k\delta_{i+ i_p}[T^k]+
+\beta_i^k\delta_{i+ i_p}[S^k]$ in terms of  $-\alpha_i^k \delta_{k+
+  k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]$.
+Where the coordinates slope, the $i$-flux gives a PE change
+\begin{multline}
+  \label{eq:triad:lat_densityPE}
+ g \delta_{i+i_p}[z_T^k]
+\left[
+-\alpha _i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (T) + \beta_i^k {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (S)
+\right] \\
+= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
+     \frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}
+\left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)
+\frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}},
+\end{multline}
+(using \eqref{eq:triad:skewfluxu}) and so the total PE change
+\eqref{eq:triad:vert_densityPE} + \eqref{eq:triad:lat_densityPE} associated with the triad fluxes is
+\begin{multline}
+  \label{eq:triad:tot_densityPE}
+  g{e_{3w}}_{\,i}^{\,k+k_p}{\mathbb{S}_w}_{i_p}^{k_p} (\rho) +
+g\delta_{i+i_p}[z_T^k] {\:}_i^k {\mathbb{S}_u}_{i_p}^{k_p} (\rho) \\
+= +\quarter g{A_e}_i^k{\: }{b_u}_{i+i_p}^k
+     \left({_i^k\mathbb{R}_{i_p}^{k_p}}+\frac{\delta_{i+i_p}[z_T^k]}{{e_{1u}}_{\,i + i_p}^{\,k}}\right)^2
+\frac{-\alpha_i^k \delta_{k+ k_p}[T^i]+ \beta_i^k\delta_{k+ k_p}[S^i]} {{e_{3w}}_{\,i}^{\,k+k_p}}.
+\end{multline}
+Where the fluid is stable, with $-\alpha_i^k \delta_{k+ k_p}[T^i]+
+\beta_i^k\delta_{k+ k_p}[S^i]<0$, this PE change is negative.
+\subsection{Treatment of the triads at the boundaries}\label{sec:triad:skew_bdry}
+Triad slopes \rtriadt{R} used for the calculation of the eddy-induced skew-fluxes
+are masked at the boundaries in exactly the same way as are the triad
+slopes \rtriad{R} used for the iso-neutral diffusive fluxes, as
+described in \S\ref{sec:triad:iso_bdry} and
+Fig.~\ref{fig:triad:bdry_triads}. Thus surface layer triads
+$\triadt{i}{1}{R}{1/2}{-1/2}$ and $\triadt{i+1}{1}{R}{-1/2}{-1/2}$ are
+masked, and both near bottom triad slopes $\triadt{i}{k}{R}{1/2}{1/2}$
+and $\triadt{i+1}{k}{R}{-1/2}{1/2}$ are masked when either of the
+$i,k+1$ or $i+1,k+1$ tracer points is masked, i.e.\ the $i,k+1$
+$u$-point is masked. The namelist parameter \np{ln\_botmix\_grif} has
+no effect on the eddy-induced skew-fluxes.
+\subsection{ Limiting of the slopes within the interior}\label{sec:triad:limitskew}
+Presently, the iso-neutral slopes $\tilde{r}_i$ relative
+to geopotentials are limited to be less than $1/100$, exactly as in
+calculating the iso-neutral diffusion, \S \ref{sec:triad:limit}. Each
+individual triad \rtriadt{R} is so limited.
+\subsection{Tapering within the surface mixed layer}\label{sec:triad:taperskew}
+The slopes $\tilde{r}_i$ relative to
+geopotentials (and thus the individual triads \rtriadt{R}) are always tapered linearly from their value immediately below the mixed layer to zero at the
+surface \eqref{eq:triad:rmtilde}, as described in \S\ref{sec:triad:lintaper}. This is
+option (c) of Fig.~\ref{Fig_eiv_slp}. This linear tapering for the
+slopes used to calculate the eddy-induced fluxes is
+unaffected by the value of \np{ln\_triad\_iso}.
+The justification for this linear slope tapering is that, for $A_e$
+that is constant or varies only in the horizontal (the most commonly
+used options in \NEMO: see \S\ref{LDF_coef}), it is
+equivalent to a horizontal eiv (eddy-induced velocity) that is uniform
+within the mixed layer \eqref{eq:triad:eiv_v}. This ensures that the
+eiv velocities do not restratify the mixed layer \citep{Treguier1997,
+  Danabasoglu_al_2008}. 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, and so a tracer convergence uniform in depth (the
+horizontal flux convergence is relatively insignificant within the mixed-layer).
+\subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag}
+Where the namelist parameter \np{ln\_traldf\_gdia}=true, diagnosed
+mean eddy-induced velocities are output. Each time step,
+streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
+$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
+\ref{Tab_cell}) respectively. We follow \citep{Griffies_Bk04} and
+calculate the streamfunction at a given $uw$-point from the
+surrounding four triads according to:
+\begin{equation}
+  \label{eq:triad:sfdiagi}
+  {\psi_1}_{i+1/2}^{k+1/2}={\quarter}\sum_{\substack{i_p,\,k_p}}
+  {A_e}_{i+1/2-i_p}^{k+1/2-k_p}\:\triadd{i+1/2-i_p}{k+1/2-k_p}{R}{i_p}{k_p}.
+\end{equation}
+The streamfunction $\psi_1$ is calculated similarly at $vw$ points.
+The eddy-induced velocities are then calculated from the
+straightforward discretisation of \eqref{eq:triad:eiv_v}:
+\begin{equation}\label{eq:triad:eiv_v_discrete}
+\begin{split}
+ {u^*}_{i+1/2}^{k} & = - \frac{1}{{e_{3u}}_{i}^{k}}\left({\psi_1}_{i+1/2}^{k+1/2}-{\psi_1}_{i+1/2}^{k+1/2}\right),   \\
+ {v^*}_{j+1/2}^{k} & = - \frac{1}{{e_{3v}}_{j}^{k}}\left({\psi_2}_{j+1/2}^{k+1/2}-{\psi_2}_{j+1/2}^{k+1/2}\right),   \\
+ {w^*}_{i,j}^{k+1/2} & =    \frac{1}{e_{1t}e_{2t}}\; \left\{
+ {e_{2u}}_{i+1/2}^{k+1/2} \,{\psi_1}_{i+1/2}^{k+1/2} -
+ {e_{2u}}_{i-1/2}^{k+1/2} \,{\psi_1}_{i-1/2}^{k+1/2} \right. + \\
+\phantom{=} & \qquad\qquad\left. {e_{2v}}_{j+1/2}^{k+1/2} \,{\psi_2}_{j+1/2}^{k+1/2} - {e_{2v}}_{j-1/2}^{k+1/2} \,{\psi_2}_{j-1/2}^{k+1/2} \right\},
+\end{split}
+\end{equation}

trunk/DOC/TexFiles/Chapters/Chap_ASM.tex

-                      r2483
+                      r3294
 The ASM code adds the functionality to apply increments to the model variables:
 temperature, salinity, sea surface height, velocity and sea ice concentration.
 These are read into the model from a NetCDF file which may be produced by data
 assimilation.  The code can also output model background fields which are used
+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. This is all controlled by the namelist
 \textit{nam\_asminc}.  There is a brief description of all the namelist options
 …
 %==========================================================================
+% Divergence damping description %%%
+\section{Divergence damping initialisation}
+\label{ASM_details}
+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:
+\begin{equation} \label{eq:asm_dmp}
+\left\{ \begin{aligned}
+ u^{n}_I = u^{n-1}_I + \frac{1}{e_{1u} } \delta _{i+1/2} \left( {A_D
+\;\chi^{n-1}_I } \right) \\
+\\
+ v^{n}_I = v^{n-1}_I + \frac{1}{e_{2v} } \delta _{j+1/2} \left( {A_D
+\;\chi^{n-1}_I } \right) \\
+\end{aligned} \right.,
+\end{equation}
+where
+\begin{equation} \label{eq:asm_div}
+\chi^{n-1}_I = \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
+                \left( {\delta _i \left[ {e_{2u}\,e_{3u}\,u^{n-1}_I} \right]
+                       +\delta _j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right).
+\end{equation}
+By the application of \eqref{eq:asm_dmp} and \eqref{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 the divergence is strongly damped.
+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}. 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.
+%==========================================================================
 \section{Implementation details}

trunk/DOC/TexFiles/Chapters/Chap_CFG.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 at $\sim 30\deg$N and rotated by 45\deg, 3180~km long, 2120~km wide
 and 4~km deep (Fig.~\ref{Fig_MISC_strait_hand}).
 The domain is bounded by vertical walls and by a ßat bottom. The configuration is
+The domain is bounded by vertical walls and by a flat bottom. 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 ßuxes.
+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}.
 …
 It forces a subpolar gyre in the north, a subtropical gyre in the wider part of the domain
 and a small recirculation gyre in the southern corner.
 The net heat ßux takes the form of a restoring toward a zonal apparent air
 temperature profile. A portion of the net heat ßux which comes from the solar radiation
+The net heat flux takes the form of a restoring toward a zonal apparent air
+temperature profile. A portion of the net heat flux which comes from the solar radiation
 is allowed to penetrate within the water column.
 The fresh water ßux is also prescribed and varies zonally.
 It is determined such as, at each time step, the basin-integrated ßux is zero.
+The fresh water flux is also prescribed and varies zonally.
+It is determined such as, at each time step, the basin-integrated flux is zero.
 The basin is initialised at rest with vertical profiles of temperature and salinity
 uniformly applied to the whole domain.
 …
 % -------------------------------------------------------------------------------------------------------------
+%       POMME configuration
+% -------------------------------------------------------------------------------------------------------------
+\section{POMME: mid-latitude sub-domain}
+\label{MISC_config_POMME}
+\key{pomme\_r025} : to be described....
+%       AMM configuration
+% -------------------------------------------------------------------------------------------------------------
+\section{AMM: atlantic margin configuration (\key{amm\_12km})}
+\label{MISC_config_AMM}
+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 key \key{amm\_12km}
+is used to create 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 (key\_bdy).
+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 addition to the tidal boundary condition the model may also take
+open boundary conditions from a North Atlantic model. Boundaries may be
+completely ommited by removing the BDY key (key\_bdy).
+Sample surface fluxes, river forcing and a sample initial restart file
+are included to test a realistic model run. The Baltic boundary is
+included within the river input file and is specified as a river source.
+Unlike ordinary river points the Baltic inputs also include salinity and
+temperature data.

trunk/DOC/TexFiles/Chapters/Chap_Conservation.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Chap_DIA.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 \item[group]: define a group of file or field. Accept the same attributes as file or field.
 \item[file]: define the output fileÕs characteristics. Accepted attributes are description, enable,
+\item[file]: define the output file's characteristics. Accepted attributes are description, enable,
 output\_freq, output\_level, id, name, name\_suffix. Child of file\_definition or group of files tag.
 …
 \item[output\_level]: file attribute. Integer from 0 to 10 defining the output priority of variables in a file:
 all variables listed in the file with a level smaller or equal to output\_level will be output.
 Other variables wonÕt be output even if they are listed in the file.
+Other variables won't be output even if they are listed in the file.
 \item[positive]: axis attribute (always .FALSE.). Logical defining the vertical axis convention used
 …
 numeric of the code, so that the trajectories never intercept the bathymetry.
+\subsubsection{ Input data: initial coordinates }
+Initial coordinates can be given with Ariane Tools convention ( IJK coordinates ,(\np{ln\_ariane}=true) )
+or with longitude and latitude.
+In case of Ariane convention, input filename is \np{"init\_float\_ariane"}. Its format is:
+\texttt{ I J K nisobfl itrash itrash }
+\noindent with:
+ - I,J,K  : indexes of initial position
+ - nisobfl: 0 for an isobar float, 1 for a float following the w velocity
+ - itrash : set to zero; it is a dummy variable to respect Ariane Tools convention
+ - itrash :set to zero; it is a dummy variable to respect Ariane Tools convention
+\noindent Example:\\
+\noindent \texttt{ 100.00000  90.00000  -1.50000 1.00000   0.00000}\\
+\texttt{ 102.00000  90.00000  -1.50000 1.00000   0.00000}\\
+\texttt{ 104.00000  90.00000  -1.50000 1.00000   0.00000}\\
+\texttt{ 106.00000  90.00000  -1.50000 1.00000   0.00000}\\
+\texttt{ 108.00000  90.00000  -1.50000 1.00000   0.00000}\\
+In the other case ( longitude and latitude ), input filename is \np{"init\_float"}. Its format is:
+\texttt{ Long Lat depth nisobfl ngrpfl itrash}
+\noindent with:
+ - Long, Lat, depth  : Longitude, latitude, depth
+ - nisobfl: 0 for an isobar float, 1 for a float following the w velocity
+ - ngrpfl : number to identify searcher group
+ - itrash :set to 1; it is a dummy variable.
+\noindent Example:
+\noindent\texttt{  20.0 0.0 0.0 0 1 1 }\\
+\texttt{ -21.0 0.0 0.0 0 1 1 }\\
+\texttt{ -22.0 0.0 0.0 0 1 1 }\\
+\texttt{ -23.0 0.0 0.0 0 1 1 }\\
+\texttt{ -24.0 0.0 0.0 0 1 1 }\\
+\np{jpnfl} is the total number of floats during the run.
+When initial positions are read in a restart file ( \np{ln\_rstflo= .TRUE.} ),  \np{jpnflnewflo}
+can be added in the initialization file.
+\subsubsection{ Output data }
+\np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl}
+is the frequency of creation of the float restart file.
+Output data can be written in ascii files (\np{ln\_flo\_ascii = .TRUE.} ). In that case,
+output filename is \np{is trajec\_float}.
+Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii = .FALSE.} ). There are 2 possibilities:
+ - if (\key{iomput}) is used, outputs are selected in  \np{iodef.xml}. Here it is an example of specification
+   to put in files description section:
+\vspace{-30pt}
+\begin{alltt}  {{\scriptsize
+\begin{verbatim}
+     <group id="1d\_grid\_T" name="auto" description="ocean T grid variables" >   }
+       <file id="floats"  description="floats variables"> }\\
+           <field ref="traj\_lon"   name="floats\_longitude"   freq\_op="86400" />}
+           <field ref="traj\_lat"   name="floats\_latitude"    freq\_op="86400" />}
+           <field ref="traj\_dep"   name="floats\_depth"       freq\_op="86400" />}
+           <field ref="traj\_temp"  name="floats\_temperature" freq\_op="86400" />}
+           <field ref="traj\_salt"  name="floats\_salinity"    freq\_op="86400" />}
+           <field ref="traj\_dens"  name="floats\_density"     freq\_op="86400" />}
+           <field ref="traj\_group" name="floats\_group"       freq\_op="86400" />}
+       </file>}
+  </group>}
+\end{verbatim}
+}}\end{alltt}
+ -  if (\key{iomput}) is not used, a file called \np{trajec\_float.nc} 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 (\key{diaharm}) }
+\label{DIA_diag_harm}
+A module is available to compute the amplitude and phase for tidal waves.
+This diagnostic is actived with \key{diaharm}.
+%------------------------------------------namdia_harm----------------------------------------------------
+\namdisplay{namdia_harm}
+%----------------------------------------------------------------------------------------------------------
+Concerning the on-line Harmonic analysis, some parameters are available in namelist:
+- \texttt{nit000\_han} is the first time step used for harmonic analysis
+- \texttt{nitend\_han} is the last time step used for harmonic analysis
+- \texttt{nstep\_han} is the time step frequency for harmonic analysis
+- \texttt{nb\_ana} is the number of harmonics to analyse
+- \texttt{tname} is an array with names of tidal constituents to analyse
+\texttt{nit000\_han} and \texttt{nitend\_han} must be between \texttt{nit000} and \texttt{nitend} of the simulation.
+The restart capability is not implemented.
+The Harmonic analysis solve this equation:
+\begin{equation}
+h_{i} - A_{0} + \sum^{nb\_ana}_{j=1}[A_{j}cos(\nu_{j}t_{j}-\phi_{j})] = e_{i}
+\end{equation}
+With $A_{j}$,$\nu_{j}$,$\phi_{j}$, the amplitude, frequency and phase for each wave and $e_{i}$ the error.
+$h_{i}$ is the sea level for the time $t_{i}$ and $A_{0}$ is the mean sea level. \\
+We can rewrite this equation:
+\begin{equation}
+h_{i} - A_{0} + \sum^{nb\_ana}_{j=1}[C_{j}cos(\nu_{j}t_{j})+S_{j}sin(\nu_{j}t_{j})] = e_{i}
+\end{equation}
+with $A_{j}=\sqrt{C^{2}_{j}+S^{2}_{j}}$ et $\phi_{j}=arctan(S_{j}/C_{j})$.
+We obtain in output $C_{j}$ and $S_{j}$ for each tidal wave.
+% -------------------------------------------------------------------------------------------------------------
+%       Sections transports
+% -------------------------------------------------------------------------------------------------------------
+\section{Transports across sections (\key{diadct}) }
+\label{DIA_diag_dct}
+A module is available to compute the transport of volume, heat and salt through sections. This diagnostic
+is actived with \key{diadct}.
+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 NEMO to compute on-line transports.
+The on-line transports module creates three output ascii files:
+- \texttt{volume\_transport} for volume transports (  unit: $10^{6} m^{3} s^{-1}$ )
+- \texttt{heat\_transport}   for heat transports   (  unit: $10^{15} W $ )
+- \texttt{salt\_transport}   for salt transports   (  unit: $10^{9}g s^{-1}$ )\\
+Namelist parameters control how frequently the flows are summed and the time scales over which
+ they are averaged, as well as the level of output for debugging:
+%------------------------------------------namdct----------------------------------------------------
+\namdisplay{namdct}
+%-------------------------------------------------------------------------------------------------------------
+\texttt{nn\_dct}: frequency of instantaneous transports computing
+\texttt{nn\_dctwri}: frequency of writing ( mean of instantaneous transports )
+\texttt{nn\_debug}: debugging of the section
+\subsubsection{ To create a binary file containing the pathway of each section }
+In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \texttt{ {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).
+Another file is available for the GYRE configuration (\texttt{ {list\_sections.ascii\_GYRE}}).
+Each section is defined by:
+\noindent \texttt{ long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name }\\
+with:
+- \texttt{long1 lat1} , coordinates of the first extremity of the section;
+- \texttt{long2 lat2} , coordinates of the second extremity of the section;
+- \texttt{nclass} the number of bounds of your classes (e.g. 3 bounds for 2 classes);
+- \texttt{okstrpond} to compute heat and salt transport, \texttt{nostrpond} if no;
+- \texttt{ice}  to compute surface and volume ice transports, \texttt{noice} if no. \\
+\noindent The results of the computing of transports, and the directions of positive
+ and negative flow do not depend on the order of the 2 extremities in this file.\\
+\noindent If nclass =/ 0,the next lines contain the class type and the nclass bounds:
+\texttt{long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name}
+\texttt{classtype}
+\texttt{zbound1}
+\texttt{zbound2}
+\texttt{.}
+\texttt{.}
+\texttt{nclass-1}
+\texttt{nclass}
+\noindent where \texttt{classtype} can be:
+- \texttt{zsal} for salinity classes
+- \texttt{ztem} for temperature classes
+- \texttt{zlay} for depth classes
+- \texttt{zsigi} for insitu density classes
+- \texttt{zsigp} for potential density classes \\
+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. \\
+It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated
+with the coordinates file name and path. \\
+Examples of two sections, the ACC\_Drake\_Passage with no classes, and the
+ ATL\_Cuba\_Florida with 4 temperature clases (5 class bounds), are shown:
+\noindent \texttt{ -68.    -54.5   -60.    -64.7  00 okstrpond noice ACC\_Drake\_Passage}
+\noindent \texttt{ -80.5    22.5   -80.5    25.5  05 nostrpond noice ATL\_Cuba\_Florida}
+\noindent \texttt{ztem}
+\noindent \texttt{-2.0}
+\noindent \texttt{ 4.5}
+\noindent \texttt{ 7.0}
+\noindent \texttt{12.0}
+\noindent \texttt{40.0}
+\subsubsection{ To read the output files }
+The output format is :
+{\small\texttt{date, time-step number, section number, section name, section slope coefficient, class number,
+class name, class bound 1 , classe bound2, transport\_direction1 ,  transport\_direction2, transport\_total}}\\
+For sections with classes, the first \texttt{nclass-1} lines correspond to the transport for each class
+and the last line corresponds to the total transport summed over all classes. For sections with no classes, class number
+\texttt{ 1 } corresponds to \texttt{ total class } and this class is called  \texttt{N}, meaning \texttt{none}.\\
+\noindent \texttt{ transport\_direction1 } is the positive part of the transport ( \texttt{ > = 0 } ).
+\noindent \texttt{ transport\_direction2 } is the negative part of the transport ( \texttt{ < = 0 } ).\\
+\noindent The \texttt{section slope coefficient} gives information about the significance of transports signs and direction:\\
+\begin{tabular}{|c|c|c|c|p{4cm}|}
+\hline
+section slope coefficient & section type & direction 1 & direction 2 & total transport \\ \hline
+.    &  horizontal & northward & southward & postive: northward  \\ \hline
+. &  vertical   & eastward  & westward  & postive: eastward  \\ \hline
+\texttt{=/0, =/ 1000.}   &  diagonal   & eastward  & westward  & postive: eastward  \\ \hline
+\end{tabular}
 % -------------------------------------------------------------------------------------------------------------
 …
 are removed from the sub-basins. Note also that the Arctic Ocean has been split
 into Atlantic and Pacific basins along the North fold line.  }
 \end{center}   \end{figure}
+\end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 (see Section \ref{MISC_steric} below for one of them).
 Activating those outputs requires to define the \key{diaar5} CPP key.
+\\
+\\

trunk/DOC/TexFiles/Chapters/Chap_DOM.tex

Property svn:executable deleted

-                      r2376
+                      r3294
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{table}[!tb] \label{Tab_cell}
+\begin{table}[!tb]
 \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
 \hline
 …
 fw & $i+1/2$   & $j+1/2$   & $k+1/2$   \\ \hline
 \end{tabular}
 \caption{ \label{Tab_cell}
+\caption{ \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.
 …
 the following discrete forms in the curvilinear $s$-coordinate system:
 \begin{equation} \label{Eq_DOM_grad}
 \nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,{\rm {\bf i}}
       +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,{\rm {\bf j}}
       +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,{\rm {\bf k}}
+\nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,\mathbf{i}
+      +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,\mathbf{j}
+      +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,\mathbf{k}
 \end{equation}
 \begin{multline} \label{Eq_DOM_lap}
 \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
+\Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
        \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right]
 +                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)      \\
 …
 \begin{eqnarray}  \label{Eq_DOM_curl}
  \nabla \times {\rm {\bf A}}\equiv &
       \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \rm{\bf i} \\
  +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \rm{\bf j} \\
  +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \rm{\bf k}
+      \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\
+ +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \mathbf{j} \\
+ +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \mathbf{k}
  \end{eqnarray}
 \begin{equation} \label{Eq_DOM_div}
 …
 \end{align*}
-\gmcomment{ STEVEN: are the dots multiplications?}
 Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with
 the numerical indexing used (\S~\ref{DOM_Num_Index}). Moreover, the
 …
 \gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  }
 %%%
+% ================================================================
+% Domain: Initial State (dtatsd & istate)
+% ================================================================
+\section  [Domain: Initial State (\textit{istate and dtatsd})]
+      {Domain: Initial State \small{(\mdl{istate} and \mdl{dtatsd} modules)} }
+\label{DTA_tsd}
+%-----------------------------------------namtsd-------------------------------------------
+\namdisplay{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.
+\begin{description}
+\item[ln\_tsd\_init = .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 \S~\ref{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[ln\_tsd\_init = .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.
+\end{description}

trunk/DOC/TexFiles/Chapters/Chap_DYN.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 the relative vorticity term and horizontal kinetic energy for the planetary vorticity
 term (MIX scheme) ; or conserving both the potential enstrophy of horizontally non-divergent
+flow and horizontal kinetic energy (ENE scheme) (see  Appendix~\ref{Apdx_C_vor_zad}).
+flow and horizontal kinetic energy (EEN scheme) (see  Appendix~\ref{Apdx_C_vor_zad}). In the
+case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the
+consistency of vorticity term with analytical equations (\textit{ln\_dynvor\_con}=true).
 The vorticity terms are all computed in dedicated routines that can be found in
 the \mdl{dynvor} module.
 …
 Pressure gradient formulations in an $s$-coordinate have been the subject of a vast
 number of papers ($e.g.$, \citet{Song1998, Shchepetkin_McWilliams_OM05}).
+A number of different pressure gradient options are coded, but they are not yet fully
+documented or tested.
+$\bullet$ Traditional coding (see for example \citet{Madec_al_JPO96}: (\np{ln\_dynhpg\_sco}=true,
+\np{ln\_dynhpg\_hel}=true)
+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}=true)
 \begin{equation} \label{Eq_dynhpg_sco}
 \left\{ \begin{aligned}
 …
 \eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of
 the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point
+($e_{3w}$). The version \np{ln\_dynhpg\_hel}=true has been added by Aike
+Beckmann and involves a redefinition of the relative position of $T$-points relative
+to $w$-points.
+$\bullet$ Weighted density Jacobian (WDJ) \citep{Song1998} (\np{ln\_dynhpg\_wdj}=true)
+($e_{3w}$).
+$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true)
 $\bullet$ Density Jacobian with cubic polynomial scheme (DJC) \citep{Shchepetkin_McWilliams_OM05}
+(\np{ln\_dynhpg\_djc}=true)
+$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true)
+Note that expression \eqref{Eq_dynhpg_sco} is 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 other pressure gradient options are not yet available.
+(\np{ln\_dynhpg\_djc}=true) (currently disabled; under development)
+Note that expression \eqref{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 pressure jacobian scheme
+(\np{ln\_dynhpg\_prj}=true) is available as an improved option to \np{ln\_dynhpg\_sco}=true when
+\key{vvl} is active.  The pressure Jacobian scheme uses a constrained cubic spline to reconstruct
+the density profile across the water column. This method maintains the monotonicity between the
+density nodes  The pressure can be calculated by analytical integration of the density profile and a
+pressure Jacobian method is used to solve the horizontal pressure gradient. This method can provide
+a more accurate calculation of the horizontal pressure gradient than the standard scheme.
 %--------------------------------------------------------------------------------------------------------------
 …
 % ================================================================
+% Neptune effect
+% ================================================================
+\section  [Neptune effect (\textit{dynnept})]
+                {Neptune effect (\mdl{dynnept})}
+\label{DYN_nept}
+The "Neptune effect" (thus named in \citep{HollowayOM86}) is a
+parameterisation of the potentially large effect of topographic form stress
+(caused by eddies) in driving the ocean circulation. Originally developed for
+low-resolution models, in which it was applied via a Laplacian (second-order)
+diffusion-like term in the momentum equation, it can also be applied in eddy
+permitting or resolving models, in which a more scale-selective bilaplacian
+(fourth-order) implementation is preferred. This mechanism has a
+significant effect on boundary currents (including undercurrents), and the
+upwelling of deep water near continental shelves.
+The theoretical basis for the method can be found in
+\citep{HollowayJPO92}, including the explanation of why form stress is not
+necessarily a drag force, but may actually drive the flow.
+\citep{HollowayJPO94} demonstrate the effects of the parameterisation in
+the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees.
+\citep{HollowayOM08} demonstrate the biharmonic version of the
+parameterisation in a global run of the POP model, with an average horizontal
+grid spacing of about 32km.
+The NEMO implementation is a simplified form of that supplied by
+Greg Holloway, the testing of which was described in \citep{HollowayJGR09}.
+The major simplification is that a time invariant Neptune velocity
+field is assumed.  This is computed only once, during start-up, and
+made available to the rest of the code via a module.  Vertical
+diffusive terms are also ignored, and the model topography itself
+is used, rather than a separate topographic dataset as in
+\citep{HollowayOM08}.  This implementation is only in the iso-level
+formulation, as is the case anyway for the bilaplacian operator.
+The velocity field is derived from a transport stream function given by:
+\begin{equation} \label{Eq_dynnept_sf}
+\psi = -fL^2H
+\end{equation}
+where $L$ is a latitude-dependant length scale given by:
+\begin{equation} \label{Eq_dynnept_ls}
+L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right )
+\end{equation}
+where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively.
+Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as:
+\begin{equation} \label{Eq_dynnept_vel}
+u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \  ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x}
+\end{equation}
+\smallskip
+%----------------------------------------------namdom----------------------------------------------------
+\namdisplay{namdyn_nept}
+%--------------------------------------------------------------------------------------------------------
+\smallskip
+The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false).
+\np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied
+to the Neptune effect flow field (default=false) (this smoothing method is as
+used by Holloway).  \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and
+polar values respectively of the length-scale parameter $L$ used in determining
+the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}.
+Values at intermediate latitudes are given by a cosine fit, mimicking the
+variation of the deformation radius with latitude.  The default values of 12km
+and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse
+resolution model. The finer resolution study of \citep{HollowayOM08} increased
+the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the
+stream function for a given topography.
+The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities
+in shallow water, and \citep{HollowayOM08} add an offset to the depth in the
+denominator to control this problem. In this implementation we offer instead (at
+the suggestion of G. Madec) the option of ramping down the Neptune flow field to
+zero over a finite depth range. The switch \np{ln\_neptramp} activates this
+option (default=false), in which case velocities at depths greater than
+\np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a
+depth of \np{rn\_htrmin} (and shallower).
+% ================================================================

trunk/DOC/TexFiles/Chapters/Chap_LBC.tex

Property svn:executable deleted

-                      r3230
+                      r3294
 %-----------------------------------------nambdy--------------------------------------------
-%-    cn_mask    =  ''                        !  name of mask file (if ln_bdy_mask=.TRUE.)
-%-    cn_dta_frs_T  = 'bdydata_grid_T.nc'     !  name of data file (T-points)
-%-    cn_dta_frs_U  = 'bdydata_grid_U.nc'     !  name of data file (U-points)
-%-    cn_dta_frs_V  = 'bdydata_grid_V.nc'     !  name of data file (V-points)
-%-    cn_dta_fla_T  = 'bdydata_bt_grid_T.nc'  !  name of data file for Flather condition (T-points)
-%-    cn_dta_fla_U  = 'bdydata_bt_grid_U.nc'  !  name of data file for Flather condition (U-points)
-%-    cn_dta_fla_V  = 'bdydata_bt_grid_V.nc'  !  name of data file for Flather condition (V-points)
-%-    ln_clim    = .false.                    !  contain 1 (T) or 12 (F) time dumps and be cyclic
-%-    ln_vol     = .true.                     !  total volume correction (see volbdy parameter)
-%-    ln_mask    = .false.                    !  boundary mask from filbdy_mask (T) or boundaries are on edges of domain (F)
-%-    ln_tides   = .true.                     !  Apply tidal harmonic forcing with Flather condition
-%-    ln_dyn_fla = .true.                     !  Apply Flather condition to velocities
-%-    ln_tra_frs = .false.                    !  Apply FRS condition to temperature and salinity
-%-    ln_dyn_frs = .false.                    !  Apply FRS condition to velocities
-%-    nn_rimwidth    =  9                     !  width of the relaxation zone
-%-    nn_dtactl      =  1                     !  = 0, bdy data are equal to the initial state
-%-                                            !  = 1, bdy data are read in 'bdydata   .nc' files
-%-    nn_volctl      =  0                     !  = 0, the total water flux across open boundaries is zero
-%-                                            !  = 1, the total volume of the system is conserved
 \namdisplay{nambdy}
+%-----------------------------------------------------------------------------------------------
+%-----------------------------------------nambdy_index--------------------------------------------
+\namdisplay{nambdy_index}
+%-----------------------------------------------------------------------------------------------
+%-----------------------------------------nambdy_dta--------------------------------------------
+\namdisplay{nambdy_dta}
+%-----------------------------------------------------------------------------------------------
+%-----------------------------------------nambdy_dta--------------------------------------------
+\namdisplay{nambdy_dta2}
 %-----------------------------------------------------------------------------------------------
 …
 The BDY module was modelled on the OBC module and shares many features
 and a similar coding structure \citep{Chanut2005}.
+The BDY module is completely rewritten at NEMO 3.4 and there is a new
+set of namelists. Boundary data files used with earlier versions of
+NEMO may need to be re-ordered to work with this version. See the
+section on the Input Boundary Data Files for details.
+%----------------------------------------------
+\subsection{The namelists}
+\label{BDY_namelist}
+It is possible to define more than one boundary ``set'' and apply
+different boundary conditions to each set. The number of boundary
+sets is defined by \np{nb\_bdy}.  Each boundary set may be defined
+as a set of straight line segments in a namelist
+(\np{ln\_coords\_file}=.false.) or read in from a file
+(\np{ln\_coords\_file}=.true.). If the set is defined in a namelist,
+then the namelists nambdy\_index must be included separately, one for
+each set. If the set is defined by a file, then a
+``coordinates.bdy.nc'' file must be provided. The coordinates.bdy file
+is analagous to the usual NEMO ``coordinates.nc'' file. In the example
+above, there are two boundary sets, the first of which is defined via
+a file and the second is defined in a namelist. For more details of
+the definition of the boundary geometry see section
+\ref{BDY_geometry}.
+For each boundary set a boundary
+condition has to be chosen for the barotropic solution (``u2d'':
+sea-surface height and barotropic velocities), for the baroclinic
+velocities (``u3d''), and for the active tracers\footnote{The BDY
+  module does not deal with passive tracers at this version}
+(``tra''). For each set of variables there is a choice of algorithm
+and a choice for the data, eg. for the active tracers the algorithm is
+set by \np{nn\_tra} and the choice of data is set by
+\np{nn\_tra\_dta}.
+The choice of algorithm is currently as follows:
+\mbox{}
+\begin{itemize}
+\item[0.] No boundary condition applied. So the solution will ``see''
+  the land points around the edge of the edge of the domain.
+\item[1.] Flow Relaxation Scheme (FRS) available for all variables.
+\item[2.] Flather radiation scheme for the barotropic variables. The
+  Flather scheme is not compatible with the filtered free surface
+  ({\it dynspg\_ts}).
+\end{itemize}
+\mbox{}
+The main choice for the boundary data is
+to use initial conditions as boundary data (\np{nn\_tra\_dta}=0) or to
+use external data from a file (\np{nn\_tra\_dta}=1). For the
+barotropic solution there is also the option to use tidal
+harmonic forcing either by itself or in addition to other external
+data.
+If external boundary data is required then the nambdy\_dta namelist
+must be defined. One nambdy\_dta namelist is required for each boundary
+set in the order in which the boundary sets are defined in nambdy. In
+the example given, two boundary sets have been defined and so there
+are two nambdy\_dta namelists. The boundary data is read in using the
+fldread module, so the nambdy\_dta namelist is in the format required
+for fldread. For each variable required, the filename, the frequency
+of the files and the frequency of the data in the files is given. Also
+whether or not time-interpolation is required and whether the data is
+climatological (time-cyclic) data. Note that on-the-fly spatial
+interpolation of boundary data is not available at this version.
+In the example namelists given, two boundary sets are defined. The
+first set is defined via a file and applies FRS conditions to
+temperature and salinity and Flather conditions to the barotropic
+variables. External data is provided in daily files (from a
+large-scale model). Tidal harmonic forcing is also used. The second
+set is defined in a namelist. FRS conditions are applied on
+temperature and salinity and climatological data is read from external
+files.
 %----------------------------------------------
 …
 Note that the sea-surface height gradient in \eqref{Eq_bdy_fla1}
 is a spatial gradient across the model boundary, so that $\eta_{e}$ is
+defined on the $T$ points with $nbrdta=1$ and $\eta$ is defined on the
+$T$ points with $nbrdta=2$. $U$ and $U_{e}$ are defined on the $U$ or
+$V$ points with $nbrdta=1$, $i.e.$ between the two $T$ grid points.
+%----------------------------------------------
+\subsection{Choice of schemes}
+\label{BDY_choice_of_schemes}
+The Flow Relaxation Scheme may be applied separately to the
+temperature and salinity (\np{ln\_tra\_frs} = true) and
+the velocity fields (\np{ln\_dyn\_frs} = true). Flather
+radiation conditions may be applied using externally defined
+barotropic velocities and sea-surface height (\np{ln\_dyn\_fla} = true)
+or using tidal harmonics fields (\np{ln\_tides} = true)
+or both. FRS and Flather conditions may be applied simultaneously.
+A typical configuration where all possible conditions might be used is a tidal,
+shelf-seas model, where the barotropic boundary conditions are fixed
+with the Flather scheme using tidal harmonics and possibly output
+from a large-scale model, and FRS conditions are applied to the tracers
+and baroclinic velocity fields, using fields from a large-scale model.
+Note that FRS conditions will work with the filtered
+(\key{dynspg\_flt}) or time-split (\key{dynspg\_ts}) solutions for the
+surface pressure gradient. The Flather condition will only work for
+the time-split solution (\key{dynspg\_ts}). FRS conditions are applied
+at the end of the main model time step. Flather conditions are applied
+during the barotropic subcycle in the time-split solution.
+defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the
+$T$ points with $nbr=2$. $U$ and $U_{e}$ are defined on the $U$ or
+$V$ points with $nbr=1$, $i.e.$ between the two $T$ grid points.
 %----------------------------------------------
 …
 \label{BDY_geometry}
+The definition of the open boundary is completely flexible. An example
+is shown in Fig.~\ref{Fig_LBC_bdy_geom}. The boundary zone is
+defined by a series of index arrays read in from the input boundary
+data files: $nbidta$, $nbjdta$, and $nbrdta$. The first two of these
+define the global $(i,j)$ indices of each point in the boundary zone
+and the $nbrdta$ array defines the discrete distance from the boundary
+with $nbrdta=1$ meaning that the point is next to the edge of the
+model domain and $nbrdta>1$ showing that the point is increasingly
+further away from the edge of the model domain. These arrays are
+defined separately for each of the $T$, $U$ and $V$ grids, but the
+relationship between the points is assumed to be as in Fig.
+\ref{Fig_LBC_bdy_geom} with the $T$ points forming the outermost row
+of the boundary and the first row of velocities normal to the boundary
+being inside the first row of $T$ points. The order in which the
+points are defined is unimportant.
+Each open boundary set is defined as a list of points. The information
+is stored in the arrays $nbi$, $nbj$, and $nbr$ in the $idx\_bdy$
+structure.  The $nbi$ and $nbj$ arrays
+define the local $(i,j)$ indices of each point in the boundary zone
+and the $nbr$ array defines the discrete distance from the boundary
+with $nbr=1$ meaning that the point is next to the edge of the
+model domain and $nbr>1$ showing that the point is increasingly
+further away from the edge of the model domain. A set of $nbi$, $nbj$,
+and $nbr$ arrays is defined for each of the $T$, $U$ and $V$
+grids. Figure \ref{Fig_LBC_bdy_geom} shows an example of an irregular
+boundary.
+The boundary geometry for each set may be defined in a namelist
+nambdy\_index or by reading in a ``coordinates.bdy.nc'' file. The
+nambdy\_index namelist defines a series of straight-line segments for
+north, east, south and west boundaries. For the northern boundary,
+\np{nbdysegn} gives the number of segments, \np{jpjnob} gives the $j$
+index for each segment and \np{jpindt} and \np{jpinft} give the start
+and end $i$ indices for each segment with similar for the other
+boundaries. These segments define a list of $T$ grid points along the
+outermost row of the boundary ($nbr\,=\, 1$). The code deduces the $U$ and
+$V$ points and also the points for $nbr\,>\, 1$ if
+$nn\_rimwidth\,>\,1$.
+The boundary geometry may also be defined from a
+``coordinates.bdy.nc'' file. Figure \ref{Fig_LBC_nc_header}
+gives an example of the header information from such a file. The file
+should contain the index arrays for each of the $T$, $U$ and $V$
+grids. The arrays must be in order of increasing $nbr$. Note that the
+$nbi$, $nbj$ values in the file are global values and are converted to
+local values in the code. Typically this file will be used to generate
+external boundary data via interpolation and so will also contain the
+latitudes and longitudes of each point as shown. However, this is not
+necessary to run the model.
+For some choices of irregular boundary the model domain may contain
+areas of ocean which are not part of the computational domain. For
+example if an open boundary is defined along an isobath, say at the
+shelf break, then the areas of ocean outside of this boundary will
+need to be masked out. This can be done by reading a mask file defined
+as \np{cn\_mask\_file} in the nam\_bdy namelist. Only one mask file is
+used even if multiple boundary sets are defined.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 \label{BDY_data}
+The input data files for the FRS conditions are defined in the
+namelist as \np{cn\_dta\_frs\_T}, \np{cn\_dta\_frs\_U},
+\np{cn\_dta\_frs\_V}. The input data files for the Flather conditions
+are defined in the namelist as \np{cn\_dta\_fla\_T},
+\np{cn\_dta\_fla\_U}, \np{cn\_dta\_fla\_V}.
+The netcdf header of a typical input data file is shown in Fig.~\ref{Fig_LBC_nc_header}.
+The file contains the index arrays which define the boundary geometry
+as noted above and the data arrays for each field.
+The data arrays are dimensioned on: a time dimension; $xb$
+which is the index of the boundary data point in the horizontal;
+and $yb$ which is a degenerate dimension of 1 to enable
+the file to be read by the standard NEMO I/O routines. The 3D fields
+also have a depth dimension.
+If \np{ln\_clim} is set to \textit{false}, the model expects the
+units of the time axis to have the form shown in
+Fig.~\ref{Fig_LBC_nc_header}, $i.e.$ {\it ``seconds since yyyy-mm-dd
+hh:mm:ss''} The fields are then linearly interpolated to the model
+time at each timestep. Note that for this option, the time axis of the
+input files must completely span the time period of the model
+integration. If \np{ln\_clim} is set to \textit{true} (climatological
+boundary forcing), the model will expect either a single set of annual
+mean fields (constant boundary forcing) or 12 sets of monthly mean
+fields in the input files.
+As in the OBC module there is an option to use initial conditions as
+boundary conditions. This is chosen by setting
+\np{nn\_dtactl}~=~0. However, since the model defines the boundary
+geometry by reading the boundary index arrays from the input files,
+it is still necessary to provide a set of input files in this
+case. They need only contain the boundary index arrays, $nbidta$,
+\textit{nbjdta}, \textit{nbrdta}.
+The data files contain the data arrays
+in the order in which the points are defined in the $nbi$ and $nbj$
+arrays. The data arrays are dimensioned on: a time dimension;
+$xb$ which is the index of the boundary data point in the horizontal;
+and $yb$ which is a degenerate dimension of 1 to enable the file to be
+read by the standard NEMO I/O routines. The 3D fields also have a
+depth dimension.
+At Version 3.4 there are new restrictions on the order in which the
+boundary points are defined (and therefore restrictions on the order
+of the data in the file). In particular:
+\mbox{}
+\begin{enumerate}
+\item The data points must be in order of increasing $nbr$, ie. all
+  the $nbr=1$ points, then all the $nbr=2$ points etc.
+\item All the data for a particular boundary set must be in the same
+  order. (Prior to 3.4 it was possible to define barotropic data in a
+  different order to the data for tracers and baroclinic velocities).
+\end{enumerate}
+\mbox{}
+These restrictions mean that data files used with previous versions of
+the model may not work with version 3.4. A fortran utility
+{\it bdy\_reorder} exists in the TOOLS directory which will re-order the
+data in old BDY data files.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_nc_header.pdf}
 \caption {     \label{Fig_LBC_nc_header}
 Example of header of netcdf input data file for BDY}
+Example of the header for a coordinates.bdy.nc file}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
 There is an option to force the total volume in the regional model to be constant,
 similar to the option in the OBC module. This is controlled  by the \np{nn\_volctl}
 parameter in the namelist. A value of\np{nn\_volctl}~=~0 indicates that this option is not used.
+parameter in the namelist. A value of \np{nn\_volctl}~=~0 indicates that this option is not used.
 If  \np{nn\_volctl}~=~1 then a correction is applied to the normal velocities
 around the boundary at each timestep to ensure that the integrated volume flow
 …
 flux across the surface and the correction velocity corrects for this as well.
+If more than one boundary set is used then volume correction is
+applied to all boundaries at once.
+\newpage
 %----------------------------------------------
 \subsection{Tidal harmonic forcing}
 \label{BDY_tides}
+%-----------------------------------------nambdy_tide--------------------------------------------
+\namdisplay{nambdy_tide}
+%-----------------------------------------------------------------------------------------------
 To be written....

trunk/DOC/TexFiles/Chapters/Chap_LDF.tex

Property svn:executable deleted

-                      r2376
+                      r3294
 and for tracers only, eddy induced advection on tracers). These three aspects
 of the lateral diffusion are set through namelist parameters and CPP keys
+(see the \textit{nam\_traldf} and \textit{nam\_dynldf} below).
+(see the \textit{nam\_traldf} and \textit{nam\_dynldf} below). Note
+that this chapter describes the default implementation of iso-neutral
+tracer mixing, and Griffies's implementation, which is used if
+\np{traldf\_grif}=true, is described in Appdx\ref{sec:triad}
 %-----------------------------------nam_traldf - nam_dynldf--------------------------------------------
 …
 $\ $\newline    % force a new ligne
 A space variation in the eddy coefficient appeals several remarks:
+The following points are relevant when the eddy coefficient varies spatially:
 (1) the momentum diffusion operator acting along model level surfaces is
 written in terms of curl and divergent components of the horizontal current
 (see \S\ref{PE_ldf}). Although the eddy coefficient can be set to different values
 in these two terms, this option is not available.
+(see \S\ref{PE_ldf}). Although the eddy coefficient could be set to different values
+in these two terms, this option is not currently available.
 (2) with an horizontally varying viscosity, the quadratic integral constraints
 …
 and either \np{ln\_traldf\_hor}=True or \np{ln\_dynldf\_hor}=True.
 \subsection{slopes for tracer iso-neutral mixing}
+\subsection{Slopes for tracer iso-neutral mixing}\label{LDF_slp_iso}
 In iso-neutral mixing  $r_1$ and $r_2$ are the slopes between the iso-neutral
 and computational surfaces. Their formulation does not depend on the vertical
 …
 \item[$s$- or hybrid $s$-$z$- coordinate : ] in the current release of \NEMO,
+there is no specific treatment for iso-neutral mixing in the $s$-coordinate.
+iso-neutral mixing is only employed for $s$-coordinates if the
+Griffies scheme is used (\np{traldf\_grif}=true; see Appdx \ref{sec:triad}).
 In other words, iso-neutral mixing will only be accurately represented with a
 linear equation of state (\np{nn\_eos}=1 or 2). In the case of a "true" equation
 …
 \end{description}
+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
+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}. Here, another strategy has been chosen \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 (Fig.~\ref{Fig_LDF_ZDF1}). This allows an iso-neutral diffusion scheme
+without additional background horizontal mixing. This technique can be viewed
+as a diffusion operator that acts along large-scale (2~$\Delta$x)
+\gmcomment{2deltax doesnt seem very large scale}
+iso-neutral surfaces. The diapycnal diffusion required for numerical stability is
+thus minimized and its net effect on the flow is quite small when compared to
+the effect of an horizontal background mixing.
+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 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 \ref{sec: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 (Fig.~\ref{Fig_LDF_ZDF1}). This allows an
+iso-neutral diffusion scheme without additional background horizontal
+mixing. This technique can be viewed as a diffusion operator that acts
+along large-scale (2~$\Delta$x) \gmcomment{2deltax doesnt seem very
+  large scale} iso-neutral surfaces. The diapycnal diffusion required
+for numerical stability is thus minimized and its net effect on the
+flow is quite small when compared to the effect of an horizontal
+background mixing.
 Nevertheless, this iso-neutral operator does not ensure that variance cannot increase,
 …
 In addition and also for numerical stability reasons \citep{Cox1987, 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).
+% In addition and also for numerical stability reasons \citep{Cox1987, 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

trunk/DOC/TexFiles/Chapters/Chap_MISC.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 Note this implementation may be sensitive to the optimization level.
+\subsection{MPP scalability}
+\label{MISC_mppsca}
+The default method of communicating values across the north-fold in distributed memory applications
+(\key{mpp\_mpi}) uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing
+region in the northern row with every other processing region in the northern row. This enables a
+global width array containing the top 4 rows to be collated on every northern row processor and then
+folded with a simple algorithm. Although conceptually simple, this "All to All" communication will
+hamper performance scalability for large numbers of northern row processors. From version 3.4
+onwards an alternative method is available which only performs direct "Peer to Peer" communications
+between each processor and its immediate "neighbours" across the fold line. This is achieved by
+using the default \textsc{MPI\_ALLGATHER} method during initialisation to help identify the "active"
+neighbours. Stored lists of these neighbours are then used in all subsequent north-fold exchanges to
+restrict exchanges to those between associated regions. The collated global width array for each
+region is thus only partially filled but is guaranteed to be set at all the locations actually
+required by each individual for the fold operation. This alternative method should give identical
+results to the default \textsc{ALLGATHER} method and is recommended for large values of \np{jpni}.
+The new method is activated by setting \np{ln\_nnogather} to be true ({\bf nammpp}). The
+reproducibility of results using the two methods should be confirmed for each new, non-reference
+configuration.
 % ================================================================

trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex

Property svn:executable deleted

-                      r2376
+                      r3294
 % Lateral Diffusive and Viscous Operators Formulation
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Lateral Diffusive and Viscous Operators Formulation}
+\subsection{Formulation of the Lateral Diffusive and Viscous Operators}
 \label{PE_ldf}
 …
 In eddy-resolving configurations, a second order operator can be used, but
 usually a more scale selective one (biharmonic operator) is preferred as the
+usually the more scale selective biharmonic operator is preferred as the
 grid-spacing is usually not small enough compared to the scale of the
 eddies. The role devoted to the subgrid-scale physics is to dissipate the
 energy that cascades toward the grid scale and thus ensures the stability of
 the model while not interfering with the solved mesoscale activity. Another approach
+energy that cascades toward the grid scale and thus to ensure the stability of
+the model while not interfering with the resolved mesoscale activity. Another approach
 is becoming more and more popular: instead of specifying explicitly a sub-grid scale
 term in the momentum and tracer time evolution equations, one uses a advective
 scheme which is diffusive enough to maintain the model stability. It must be emphasised
 that then, all the sub-grid scale physics is in this case include in the formulation of the
+that then, all the sub-grid scale physics is included in the formulation of the
 advection scheme.
 All these parameterisations of subgrid scale physics present advantages and
+All these parameterisations of subgrid scale physics have advantages and
 drawbacks. There are not all available in \NEMO. In the $z$-coordinate
 formulation, five options are offered for active tracers (temperature and
 …
 operator acting along $s-$surfaces (see \S\ref{LDF}).
 \subsubsection{lateral second order tracer diffusive operator}
+\subsubsection{Lateral second order tracer diffusive operator}
 The lateral second order tracer diffusive operator is defined by (see Appendix~\ref{Apdx_B}):
 …
 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral
+and computational surfaces. Therefore, they have a same expression in $z$- and $s$-coordinates:
+and computational surfaces. Therefore, they are different quantities,
+but have similar expressions in $z$- and $s$-coordinates. In $z$-coordinates:
 \begin{equation} \label{Eq_PE_iso_slopes}
 r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right)
                   \left( {\frac{\partial \rho }{\partial k}} \right)^{-1} \ , \quad
 r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right)
+                  \left( {\frac{\partial \rho }{\partial k}} \right)^{-1}
+\end{equation}
+When the \textit{Eddy Induced Velocity} parametrisation (eiv) \citep{Gent1990} is used,
+                  \left( {\frac{\partial \rho }{\partial k}} \right)^{-1},
+\end{equation}
+while in $s$-coordinates $\partial/\partial k$ is replaced by
+$\partial/\partial s$.
+\subsubsection{Eddy induced velocity}
+ When the \textit{eddy induced velocity} parametrisation (eiv) \citep{Gent1990} is used,
 an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers:
 \begin{equation} \label{Eq_PE_iso+eiv}
 …
 where $A^{eiv}$ is the eddy induced velocity coefficient (or equivalently the isoneutral
 thickness diffusivity coefficient), and $\tilde{r}_1$ and $\tilde{r}_2$ are the slopes
 between isoneutral and \emph{geopotential} surfaces and thus depends on the coordinate
 considered:
+between isoneutral and \emph{geopotential} surfaces. Their values are
+thus independent of the vertical coordinate, but their expression depends on the coordinate:
 \begin{align} \label{Eq_PE_slopes_eiv}
 \tilde{r}_n = \begin{cases}
 …
 to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}).
 \subsubsection{lateral fourth order tracer diffusive operator}
+\subsubsection{Lateral fourth order tracer diffusive operator}
 The lateral fourth order tracer diffusive operator is defined by:
 …
 \subsubsection{lateral second order momentum diffusive operator}
+\subsubsection{Lateral second order momentum diffusive operator}
 The second order momentum diffusive operator along $z$- or $s$-surfaces is found by

trunk/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Chap_OBS.tex

-                      r2483
+                      r3294
 $\ $\newline    % force a new line
 The observation and model comparison code (OBS) reads in observation files
 (profile temperature and salinity, sea surface temperature, sea level anomaly,
+sea ice concentration, and velocity) and calculates  an interpolated model equivalent
+value at the observation location and nearest model timestep. The OBS code is
+called from \np{opa.F90} in order to initialise the model and to calculate the
+model equivalent values for observations on the 0th timestep. The code is then
+called again after each timestep from \np{step.F90}. The code was originally
+developed for use with NEMOVAR.
+For all data types a 2D horizontal  interpolator is needed
 to interpolate the model fields to the observation location.
+For {\em in situ} profiles, a 1D vertical interpolator is needed in addition to
+provide model fields at the observation depths. Currently this only works in
+z-level model configurations, but is being developed to work with a
+generalised vertical coordinate system.
+Temperature data from moored buoys (TAO, TRITON, PIRATA) in the
+ENACT/ENSEMBLES data-base are available as daily averaged quantities. For this
+type of observation the
 observation operator will compare such observations to the model temperature
+The observation and model comparison code (OBS) reads in observation files (profile
+temperature and salinity, sea surface temperature, sea level anomaly, sea ice concentration,
+and velocity) and calculates  an interpolated model equivalent value at the observation
+location and nearest model timestep. The resulting data are saved in a ``feedback'' file (or
+files). The code was originally developed for use with the NEMOVAR data assimilation code, but
+can be used for validation or verification of model or  any other data assimilation system.
+The OBS code is called from \np{opa.F90} for model initialisation and to calculate the model
+equivalent values for observations on the 0th timestep. The code is then called again after
+each timestep from \np{step.F90}. To build with the OBS code active \key{diaobs} must be
+set.
+For all data types a 2D horizontal  interpolator is needed to interpolate the model fields to
+the observation location. For {\em in situ} profiles, a 1D vertical interpolator is needed in
+addition to provide model fields at the observation depths. Currently this only works in
+z-level model configurations, but is being developed to work with a generalised vertical
+coordinate system. Temperature data from moored buoys (TAO, TRITON, PIRATA) in the
+ENACT/ENSEMBLES data-base are available as daily averaged quantities. For this type of
+observation the observation operator will compare such observations to the model temperature
 fields averaged over one day. The relevant observation type may be specified in the namelist
+using \np{endailyavtypes}. Otherwise the model value from the nearest
+timestep to the observation time is used.
+The resulting data are saved in a ``feedback'' file (or files) which can be used
+for model validation and verification and also to provide information for data
+assimilation. This code is controlled by the namelist \textit{nam\_obs}. To
+build with the OBS code active \key{diaobs} must be set.
+Section~\ref{OBS_example} introduces a test example of the observation operator
+code including where to obtain data and how to setup the namelist.
+Section~\ref{OBS_details} introduces some more technical details of the
+different observation types used and also shows a more complete namelist.
+Finally section~\ref{OBS_theory} introduces some of the theoretical aspects of
+the observation operator including interpolation methods and running on multiple
+processors.
+using \np{endailyavtypes}. Otherwise the model value from the nearest timestep to the
+observation time is used.
+The code is controlled by the namelist \textit{nam\_obs}. See the following sections for more
+details on setting up the namelist.
+Section~\ref{OBS_example} introduces a test example of the observation operator code including
+where to obtain data and how to setup the namelist. Section~\ref{OBS_details} introduces some
+more technical details of the different observation types used and also shows a more complete
+namelist. Section~\ref{OBS_theory} introduces some of the theoretical aspects of the
+observation operator including interpolation methods and running on multiple processors.
+Section~\ref{OBS_obsutils} introduces some utilities to help working with the files produced
+by the OBS code.
 % ================================================================
 …
 \item Add the following to the NEMO namelist to run the observation
 operator on this data. Set the \np{enactfiles} namelist parameter to the
 observation  file name (or link in to \np{profiles\_01\.nc}):
+observation  file name:
 \end{enumerate}
 …
 %-------------------------------------------------------------------------------------------------------------
 The option \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
+The options \np{ln\_t3d} and \np{ln\_s3d} switch on the temperature and salinity
 profile observation operator code. The \np{ln\_ena} switch turns on the reading
 of ENACT/ENSEMBLES type profile data. The filename or array of filenames are
 …
 Setting \np{ln\_grid\_global} means that the code distributes the observations
 evenly between processors. Alternatively each processor will work with
 observations located within the model subdomain.
 The NEMOVAR system contains utilities to plot the feedback files, convert and
 recombine the files. These are available on request from the NEMOVAR team.
+observations located within the model subdomain (see section~\ref{OBS_parallel}).
+A number of utilities are now provided to plot the feedback files, convert and
+recombine the files. These are explained in more detail in section~\ref{OBS_obsutils}.
 \section{Technical details}
 …
 \subsection{Parallel aspects of horizontal interpolation}
+\label{OBS_parallel}
 For horizontal interpolation, there is the basic problem that the
 …
 \subsection{Vertical interpolation operator}
 The vertical interpolation is achieved using either a cubic spline or
+Vertical interpolation is achieved using either a cubic spline or
 linear interpolation. For the cubic spline, the top and
 bottom boundary conditions for the second derivative of the
 interpolating polynomial in the spline are set to zero.
 At the bottom boundary, this is done using the land-ocean mask.
+\newpage
+\section{Observation Utilities}
+\label{OBS_obsutils}
+Some tools for viewing and processing of observation and feedback files are provided in the
+NEMO repository for convenience. These include OBSTOOLS which are a collection of Fortran
+programs which are helpful to deal with feedback files. They do such tasks as observation file
+conversion, printing of file contents, some basic statistical analysis of feedback files. The
+other tool is an IDL program called dataplot which uses a graphical interface to visualise
+observations and feedback files. OBSTOOLS and dataplot are described in more detail below.
+\subsection{Obstools}
+A series of Fortran utilities is provided with NEMO called OBSTOOLS. This are helpful in
+handling observation files and the feedback file output from the NEMO observation operator.
+The utilities are as follows
+\subsubsection{corio2fb}
+The program corio2fb converts profile observation files from the Coriolis format to the
+standard feedback format. The program is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+corio2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+\subsubsection{enact2fb}
+The program enact2fb converts profile observation files from the ENACT format to the standard
+feedback format. The program is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+enact2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbcomb}
+The program fbcomb combines multiple feedback files produced by individual processors in an
+MPI run of NEMO into a single feedback file. The program is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbcomb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbmatchup}
+The program fbmatchup will match observations from two feedback files. The program is called
+in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbmatchup.exe outputfile inputfile1 varname1 inputfile2 varname2 ...
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbprint}
+The program fbprint will print the contents of a feedback file or files to standard output.
+Selected information can be output using optional arguments. The program is called in the
+following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbprint.exe [options] inputfile
+options:
+     -b            shorter output
+     -q            Select observations based on QC flags
+     -Q            Select observations based on QC flags
+     -B            Select observations based on QC flags
+     -u            unsorted
+     -s ID         select station ID
+     -t TYPE       select observation type
+     -v NUM1-NUM2  select variable range to print by number
+                      (default all)
+     -a NUM1-NUM2  select additional variable range to print by number
+                      (default all)
+     -e NUM1-NUM2  select extra variable range to print by number
+                      (default all)
+     -d            output date range
+     -D            print depths
+     -z            use zipped files
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbsel}
+The program fbsel will select or subsample observations. The program is called in the
+following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbsel.exe <input filename> <output filename>
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbstat}
+The program fbstat will output summary statistics in different global areas into a number of
+files. The program is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbstat.exe [-nmlev] <filenames>
+\end{verbatim}
+\end{alltt}
+\subsubsection{fbthin}
+The program fbthin will thin the data to 1 degree resolution. The code could easily be
+modified to thin to a different resolution. The program is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+fbthin.exe inputfile outputfile
+\end{verbatim}
+\end{alltt}
+\subsubsection{sla2fb}
+The program sla2fb will convert an AVISO SLA format file to feedback format. The program is
+called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+sla2fb.exe [-s type] outputfile inputfile1 inputfile2 ...
+Option:
+     -s            Select altimeter data_source
+\end{verbatim}
+\end{alltt}
+\subsubsection{vel2fb}
+The program vel2fb will convert TAO/PIRATA/RAMA currents files to feedback format. The program
+is called in the following way:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+vel2fb.exe outputfile inputfile1 inputfile2 ...
+\end{verbatim}
+\end{alltt}
+\subsection{building the obstools}
+To build the obstools use in the tools directory use ./maketools -n OBSTOOLS -m [ARCH].
+\subsection{Dataplot}
+An IDL program called dataplot is included which uses a graphical interface to visualise
+observations and feedback files. It is possible to zoom in, plot individual profiles and
+calculate some basic statistics. To plot some data run IDL and then:
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+IDL> dataplot, "filename"
+\end{verbatim}
+\end{alltt}
+To read multiple files into dataplot, for example multiple feedback files from different
+processors or from different days, the easiest method is to use the spawn command to generate
+a list of files which can then be passed to dataplot.
+\begin{alltt}
+\footnotesize
+\begin{verbatim}
+IDL> spawn, 'ls profb*.nc', files
+IDL> dataplot, files
+\end{verbatim}
+\end{alltt}
+Fig~\ref{fig:obsdataplotmain} shows the main window which is launched when dataplot starts.
+This is split into three parts. At the top there is a menu bar which contains a variety of
+drop down menus. Areas - zooms into prespecified regions; plot - plots the data as a
+timeseries or a T-S diagram if appropriate; Find - allows data to be searched; Config - sets
+various configuration options.
+The middle part is a plot of the geographical location of the observations. This will plot the
+observation value, the model background value or observation minus background value depending
+on the option selected in the radio button at the bottom of the window. The plotting colour
+range can be changed by clicking on the colour bar. The title of the plot gives some basic
+information about the date range and depth range shown, the extreme values, and the mean and
+rms values. It is possible to zoom in using a drag-box. You may also zoom in or out using the
+mouse wheel.
+The bottom part of the window controls what is visible in the plot above. There are two bars
+which select the level range plotted (for profile data). The other bars below select the date
+range shown. The bottom of the figure allows the option to plot the mean, root mean square,
+standard deviation or mean square values. As mentioned above you can choose to plot the
+observation value, the model background value or observation minus background value. The next
+group of radio buttons selects the map projection. This can either be regular latitude
+longitude grid, or north or south polar stereographic. The next group of radio buttons will
+plot bad observations, switch to salinity and plot density for profile observations. The
+rightmost group of buttons will print the plot window as a postscript, save it as png, or exit
+from dataplot.
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}     \begin{center}
+%\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_main}
+\includegraphics[width=9cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_main}
+\caption{      \label{fig:obsdataplotmain}
+Main window of dataplot.}
+\end{center}     \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+If a profile point is clicked with the mouse button a plot of the observation and background
+values as a function of depth (Fig~\ref{fig:obsdataplotprofile}).
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}     \begin{center}
+%\includegraphics[width=10cm,height=12cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_prof}
+\includegraphics[width=7cm,angle=-90.]{./TexFiles/Figures/Fig_OBS_dataplot_prof}
+\caption{      \label{fig:obsdataplotprofile}
+Profile plot from dataplot produced by right clicking on a point in the main window.}
+\end{center}     \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>

trunk/DOC/TexFiles/Chapters/Chap_SBC.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 \end{itemize}
 Four different ways to provide the first six fields to the ocean are available which
+Five different ways to provide the first six fields to the ocean are available which
 are controlled by namelist variables: an analytical formulation (\np{ln\_ana}~=~true),
 a flux formulation (\np{ln\_flx}~=~true), a bulk formulae formulation (CORE
+(\np{ln\_core}~=~true) or CLIO (\np{ln\_clio}~=~true) bulk formulae) and a coupled
+(\np{ln\_core}~=~true), CLIO (\np{ln\_clio}~=~true) or MFS
+\footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
+(\np{ln\_mfs}~=~true) bulk formulae) and a coupled
 formulation (exchanges with a atmospheric model via the OASIS coupler)
 (\np{ln\_cpl}~=~true). When used, the atmospheric pressure forces both
+ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true)
+\footnote{The surface pressure field could be use in bulk formulae, nevertheless
+none of the current bulk formulea (CLIO and CORE) uses the it.}.
+ocean and ice dynamics (\np{ln\_apr\_dyn}~=~true).
 The frequency at which the six or seven fields have to be updated is the \np{nn\_fsbc}
 namelist parameter.
 …
 (\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater
 fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of a freshwater flux adjustment
 in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); and the
+in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the
 transformation of the solar radiation (if provided as daily mean) into a diurnal
+cycle (\np{ln\_dm2dc}~=~true).
+cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave
+model (\np{ln\_cdgw}~=~true). The latter option is possible only in case core or mfs bulk formulas are selected.
 In this chapter, we first discuss where the surface boundary condition appears in the
 model equations. Then we present the four ways of providing the surface boundary condition,
+model equations. Then we present the five ways of providing the surface boundary condition,
 followed by the description of the atmospheric pressure and the river runoff.
 Next the scheme for interpolation on the fly is described.
 …
 actual year/month/day, always coded with 4 or 2 digits. Note that (1) in mpp, if the file is split
 over each subdomain, the suffix '.nc' is replaced by '\_PPPP.nc', where 'PPPP' is the
+process number coded with 4 digits; (2) when using AGRIF, the prefix ÔN\_Õ is added to files,
+process number coded with 4 digits; (2) when using AGRIF, the prefix
+'\_N' is added to files,
 where 'N'  is the child grid number.}
 \end{table}
 …
 % Bulk formulation
 % ================================================================
 \section  [Bulk formulation (\textit{sbcblk\_core} or \textit{sbcblk\_clio}) ]
       {Bulk formulation \small{(\mdl{sbcblk\_core} or \mdl{sbcblk\_clio} module)} }
+\section  [Bulk formulation (\textit{sbcblk\_core}, \textit{sbcblk\_clio} or \textit{sbcblk\_mfs}) ]
+      {Bulk formulation \small{(\mdl{sbcblk\_core} \mdl{sbcblk\_clio} \mdl{sbcblk\_mfs} modules)} }
 \label{SBC_blk}
 …
 using bulk formulae and atmospheric fields and ocean (and ice) variables.
 The atmospheric fields used depend on the bulk formulae used. Two bulk formulations
 are available : the CORE and CLIO bulk formulea. The choice is made by setting to true
 one of the following namelist variable : \np{ln\_core} and \np{ln\_clio}.
 Note : in forced mode, when a sea-ice model is used, a bulk formulation have to be used.
 Therefore the two bulk formulea provided include the computation of the fluxes over both
+The atmospheric fields used depend on the bulk formulae used. Three bulk formulations
+are available : the CORE, the CLIO and the MFS bulk formulea. The choice is made by setting to true
+one of the following namelist variable : \np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
+Note : in forced mode, when a sea-ice model is used, a bulk formulation (CLIO or CORE) have to be used.
+Therefore the two bulk (CLIO and CORE) formulea include the computation of the fluxes over both
 an ocean and an ice surface.
 …
 namelist (see \S\ref{SBC_fldread}).
+% -------------------------------------------------------------------------------------------------------------
+%        MFS Bulk formulae
+% -------------------------------------------------------------------------------------------------------------
+\subsection    [MFS Bulk formulea (\np{ln\_mfs}=true)]
+            {MFS Bulk formulea (\np{ln\_mfs}=true, \mdl{sbcblk\_mfs})}
+\label{SBC_blk_mfs}
+%------------------------------------------namsbc_mfs----------------------------------------------------
+\namdisplay{namsbc_mfs}
+%----------------------------------------------------------------------------------------------------------
+The MFS (Mediterranean Forecasting System) bulk formulae have been developed by
+ \citet{Castellari_al_JMS1998}.
+They have been designed to handle the ECMWF operational data and are currently
+in use in the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}.
+The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}.
+The surface boundary condition for temperature involves the balance between surface solar radiation,
+net long-wave radiation, the latent and sensible heat fluxes.
+Solar radiation is dependent on cloud cover and is computed by means of
+an astronomical formula \citep{Reed_JPO77}. Albedo monthly values are from \citet{Payne_JAS72}
+as means of the values at $40^{o}N$ and $30^{o}N$ for the Atlantic Ocean (hence the same latitudinal
+band of the Mediterranean Sea). The net long-wave radiation flux
+\citep{Bignami_al_JGR95} is a function of
+air temperature, sea-surface temperature, cloud cover and relative humidity.
+Sensible heat and latent heat fluxes are computed by classical
+bulk formulae parameterized according to \citet{Kondo1975}.
+Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}.
+The required 7 input fields must be provided on the model Grid-T and  are:
+\begin{itemize}
+\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
+\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
+\item          Total Claud Cover (\%)  (\np{sn\_clc})
+\item          2m Air Temperature ($K$) (\np{sn\_tair})
+\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
+\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
+\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
+\end{itemize}
+% -------------------------------------------------------------------------------------------------------------
 % ================================================================
 % Coupled formulation
 …
 \footnote{The \key{oasis4} exist. It activates portion of the code that are still under development.}.
 It has been successfully used to interface \NEMO to most of the European atmospheric
 GCM (ARPEGE, ECHAM, ECMWF, HadAM, LMDz),
+GCM (ARPEGE, ECHAM, ECMWF, HadAM, HadGAM, LMDz),
 as well as to \href{http://wrf-model.org/}{WRF} (Weather Research and Forecasting Model).
 …
 When PISCES biogeochemical model (\key{top} and \key{pisces}) is also used in the coupled system,
 the whole carbon cycle is computed by defining \key{cpl\_carbon\_cycle}. In this case,
+CO$_2$ fluxes are exchanged between the atmosphere and the ice-ocean system.
+CO$_2$ fluxes will be exchanged between the atmosphere and the ice-ocean system (and need to be activated
+in namsbc{\_}cpl).
+The new namelist above allows control of various aspects of the coupling fields (particularly for
+vectors) and now allows for any coupling fields to have multiple sea ice categories (as required by LIM3
+and CICE).  When indicating a multi-category coupling field in namsbc{\_}cpl the number of categories will be
+determined by the number used in the sea ice model.  In some limited cases it may be possible to specify
+single category coupling fields even when the sea ice model is running with multiple categories - in this
+case the user should examine the code to be sure the assumptions made are satisfactory.  In cases where
+this is definitely not possible the model should abort with an error message.  The new code has been tested using
+ECHAM with LIM2, and HadGAM3 with CICE but although it will compile with \key{lim3} additional minor code changes
+may be required to run using LIM3.
 …
 % ================================================================
+%        Tidal Potential
+% ================================================================
+\section   [Tidal Potential (\textit{sbctide})]
+                        {Tidal Potential (\mdl{sbctide})}
+\label{SBC_tide}
+A module is available to use the tidal potential forcing and is activated with with \key{tide}.
+%------------------------------------------nam_tide----------------------------------------------------
+\namdisplay{nam_tide}
+%-------------------------------------------------------------------------------------------------------------
+Concerning the tidal potential, some parameters are available in namelist:
+- \texttt{ln\_tide\_pot} activate the tidal potential forcing
+- \texttt{nb\_harmo} is the number of constituent used
+- \texttt{clname} is the name of constituent
+The tide is generated by the forces of gravity ot the Earth-Moon and Earth-Sun sytem;
+they are expressed as the gradient of the astronomical potential ($\vec{\nabla}\Pi_{a}$). \\
+The potential astronomical expressed, for the three types of tidal frequencies
+following, by : \\
+Tide long period :
+\begin{equation}
+\Pi_{a}=gA_{k}(\frac{1}{2}-\frac{3}{2}sin^{2}\phi)cos(\omega_{k}t+V_{0k})
+\end{equation}
+diurnal Tide :
+\begin{equation}
+\Pi_{a}=gA_{k}(sin 2\phi)cos(\omega_{k}t+\lambda+V_{0k})
+\end{equation}
+Semi-diurnal tide:
+\begin{equation}
+\Pi_{a}=gA_{k}(cos^{2}\phi)cos(\omega_{k}t+2\lambda+V_{0k})
+\end{equation}
+$A_{k}$ is the amplitude of the wave k, $\omega_{k}$ the pulsation of the wave k, $V_{0k}$ the astronomical phase of the wave
+$k$ to Greenwich.
+We make corrections to the astronomical potential.
+We obtain :
+\begin{equation}
+\Pi-g\delta = (1+k-h) \Pi_{A}(\lambda,\phi)
+\end{equation}
+with $k$ a number of Love estimated to 0.6 which parametrized the astronomical tidal land,
+and $h$ a number of Love to 0.3 which parametrized the parametrization due to the astronomical tidal land.
+% ================================================================
 %        River runoffs
 % ================================================================
 …
 %To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the tra_sbc module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below:
+}
+%}
 % ================================================================
 …
 ice-ocean fluxes, that are combined with the air-sea fluxes using the ice fraction of
 each model cell to provide the surface ocean fluxes. Note that the activation of a
 sea-ice model is is done by defining a CPP key (\key{lim2} or \key{lim3}).
 The activation automatically ovewrite the read value of nn{\_}ice to its appropriate
 value ($i.e.$ $2$ for LIM-2 and $3$ for LIM-3).
+sea-ice model is is done by defining a CPP key (\key{lim2}, \key{lim3} or \key{cice}).
+The activation automatically overwrites the read value of nn{\_}ice to its appropriate
+value ($i.e.$ $2$ for LIM-2, $3$ for LIM-3 or $4$ for CICE).
 \end{description}
 % {Description of Ice-ocean interface to be added here or in LIM 2 and 3 doc ?}
+\subsection   [Interface to CICE (\textit{sbcice\_cice})]
+         {Interface to CICE (\mdl{sbcice\_cice})}
+\label{SBC_cice}
+It is now possible to couple a global NEMO configuration (without AGRIF) to the CICE sea-ice
+model by using \key{cice}.  The CICE code can be obtained from
+\href{http://oceans11.lanl.gov/trac/CICE/}{LANL} and the additional 'hadgem3' drivers will be required,
+even with the latest code release.  Input grid files consistent with those used in NEMO will also be needed,
+and CICE CPP keys \textbf{ORCA\_GRID}, \textbf{CICE\_IN\_NEMO} and \textbf{coupled} should be used (seek advice from UKMO
+if necessary).  Currently the code is only designed to work when using the CORE forcing option for NEMO (with
+\textit{calc\_strair~=~true} and \textit{calc\_Tsfc~=~true} in the CICE name-list), or alternatively when NEMO
+is coupled to the HadGAM3 atmosphere model (with \textit{calc\_strair~=~false} and \textit{calc\_Tsfc~=~false}).
+The code is intended to be used with \np{nn\_fsbc} set to 1 (although coupling ocean and ice less frequently
+should work, it is possible the calculation of some of the ocean-ice fluxes needs to be modified slightly - the
+user should check that results are not significantly different to the standard case).
+There are two options for the technical coupling between NEMO and CICE.  The standard version allows
+complete flexibility for the domain decompositions in the individual models, but this is at the expense of global
+gather and scatter operations in the coupling which become very expensive on larger numbers of processors. The
+alternative option (using \key{nemocice\_decomp} for both NEMO and CICE) ensures that the domain decomposition is
+identical in both models (provided domain parameters are set appropriately, and
+\textit{processor\_shape~=~square-ice} and \textit{distribution\_wght~=~block} in the CICE name-list) and allows
+much more efficient direct coupling on individual processors.  This solution scales much better although it is at
+the expense of having more idle CICE processors in areas where there is no sea ice.
 % -------------------------------------------------------------------------------------------------------------
 …
 \end{description}
+% -------------------------------------------------------------------------------------------------------------
+%        Neutral Drag Coefficient from external wave model
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Neutral drag coefficient from external wave model (\textit{sbcwave})]
+                        {Neutral drag coefficient from external wave model (\mdl{sbcwave})}
+\label{SBC_wave}
+%------------------------------------------namwave----------------------------------------------------
+\namdisplay{namsbc_wave}
+%-------------------------------------------------------------------------------------------------------------
+\begin{description}
+\item [??] In order to read a neutral drag coeff, from an external data source (i.e. a wave model), the
+logical variable \np{ln\_cdgw}
+ in $namsbc$ namelist must be defined ${.true.}$.
+The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the
+namelist ${namsbc\_wave}$ (for external data names, locations, frequency, interpolation and all
+the miscellanous options allowed by Input Data generic Interface see \S\ref{SBC_input})
+and a 2D field of neutral drag coefficient. Then using the routine
+TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, the drag coefficient is computed according
+to stable/unstable conditions of the air-sea interface following \citet{Large_Yeager_Rep04}.
+\end{description}
 % Griffies doc:
 % When running ocean-ice simulations, we are not explicitly representing land processes, such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, it is important to balance the hydrological cycle in ocean-ice models. We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. The result of the normalization should be a global integrated zero net water input to the ocean-ice system over a chosen time scale.
 …

trunk/DOC/TexFiles/Chapters/Chap_STP.tex

Property svn:executable deleted

trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 \label{TRA_ldf_iso}
+The general form of the second order lateral tracer subgrid scale physics
+If the Griffies trad scheme is not employed
+(\np{ln\_traldf\_grif}=true; see App.\ref{sec:triad}) the general form of the second order lateral tracer subgrid scale physics
 (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and
 $s$-coordinates:
 …
 (see \S\ref{LDF}) allows the model to run safely without any additional
 background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme
 developed by \cite{Griffies_al_JPO98} which preserves both tracer and its variance
+developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
 is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of
 the algorithm is given in App.\ref{Apdx_Griffies}.
+the algorithm is given in App.\ref{sec:triad}.
 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal
 …
 %--------------------------------------------namtra_dmp-------------------------------------------------
 \namdisplay{namtra_dmp}
-\namdisplay{namdta_tem}
-\namdisplay{namdta_sal}
 %--------------------------------------------------------------------------------------------------------------
 …
 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$
 are given temperature and salinity fields (usually a climatology).
 The restoring term is added when \key{tradmp} is defined.
 It also requires that both \key{dtatem} and \key{dtasal} are defined
 and fill in \textit{namdta\_tem} and \textit{namdta\_sal namelists}
 ($i.e.$ that $T_o$ and $S_o$ are read using \mdl{fldread},
 see \S\ref{SBC_fldread}). The restoring coefficient $\gamma$ is
 a three-dimensional array initialized by the user in routine \rou{dtacof}
 also located in module \mdl{tradmp}.
+The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true.
+It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true
+in \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are
+correctly set  ($i.e.$ that $T_o$ and $S_o$ are provided in input files and read
+using \mdl{fldread}, see \S\ref{SBC_fldread}).
+The restoring coefficient $\gamma$ is a three-dimensional array initialized by the
+user in routine \rou{dtacof} also located in module \mdl{tradmp}.
 The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)}

trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex

Property svn:executable deleted

-                      r2541
+                      r3294
 % ================================================================
 % Chapter Ñ Vertical Ocean Physics (ZDF)
+% Chapter  Vertical Ocean Physics (ZDF)
 % ================================================================
 \chapter{Vertical Ocean Physics (ZDF)}
 …
 $a=5$ and $n=2$. The last three values can be modified by setting the
 \np{rn\_avmri}, \np{rn\_alp} and \np{nn\_ric} namelist parameters, respectively.
+A simple mixing-layer model to transfer and dissipate the atmospheric
+ forcings (wind-stress and buoyancy fluxes) can be activated setting
+the \np{ln\_mldw} =.true. in the namelist.
+In this case, the local depth of turbulent wind-mixing or "Ekman depth"
+ $h_{e}(x,y,t)$ is evaluated and the vertical eddy coefficients prescribed within this layer.
+This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:
+\begin{equation}
+         h_{e} = Ek \frac {u^{*}} {f_{0}}    \\
+\end{equation}
+where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis
+parameter.
+In this similarity height relationship, the turbulent friction velocity:
+\begin{equation}
+         u^{*} = \sqrt \frac {|\tau|} {\rho_o}     \\
+\end{equation}
+is computed from the wind stress vector $|\tau|$ and the reference dendity $ \rho_o$.
+The final $h_{e}$ is further constrained by the adjustable bounds \np{rn\_mldmin} and \np{rn\_mldmax}.
+Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to
+the empirical values \np{rn\_wtmix} and \np{rn\_wvmix} \citep{Lermusiaux2001}.
 % -------------------------------------------------------------------------------------------------------------
 …
 the clipping factor is of crucial importance for the entrainment depth predicted in
 stably stratified situations, and that its value has to be chosen in accordance
 with the algebraic model for the turbulent ßuxes. The clipping is only activated
+with the algebraic model for the turbulent fluxes. The clipping is only activated
 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value.
 …
 reduced as necessary to ensure stability; these changes are not reported.
+Limits on the bottom friction coefficient are not imposed if the user has elected to
+handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential
+breaches of the explicit stability criterion are still reported for information purposes.
+% -------------------------------------------------------------------------------------------------------------
+%       Implicit Bottom Friction
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Implicit Bottom Friction (\np{ln\_bfrimp}$=$\textit{T})}
+\label{ZDF_bfr_imp}
+An optional implicit form of bottom friction has been implemented to improve
+model stability. We recommend this option for shelf sea and coastal ocean applications, especially
+for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp}
+to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false}
+in the \textit{namzdf} namelist.
+This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the
+bottom boundary condition is implemented implicitly.
+\begin{equation} \label{Eq_dynzdf_bfr}
+\left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk}
+    = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}}
+\end{equation}
+where $mbk$ is the layer number of the bottom wet layer. superscript $n+1$ means the velocity used in the
+friction formula is to be calculated, so, it is implicit.
+If split-explicit time splitting is used, care must be taken to avoid the double counting of
+the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic
+pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove
+the bottom friction induced by these two terms which has been included in the 3-D momentum trend
+and update it with the latest value. On the other hand, the bottom friction contributed by the
+other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations
+and should not be added in the 2-D barotropic mode.
+The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the
+following:
+\begin{equation} \label{Eq_dynspg_ts_bfr1}
+\frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b}
+\left(\textbf{U}_{med}-\textbf{U}^{m-1}\right)
+\end{equation}
+\begin{equation} \label{Eq_dynspg_ts_bfr2}
+\frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+
+\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)-
+\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right)
+\end{equation}
+where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping
+is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step.
+ $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops
+while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom
+layer horizontal velocity.
 % -------------------------------------------------------------------------------------------------------------
 %       Bottom Friction with split-explicit time splitting
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Bottom Friction with split-explicit time splitting}
+\subsection{Bottom Friction with split-explicit time splitting (\np{ln\_bfrimp}$=$\textit{F})}
 \label{ZDF_bfr_ts}
 …
 {\key{dynspg\_flt}). Extra attention is required, however, when using
 split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface
 equation is solved with a small time step \np{nn\_baro}*\np{rn\_rdt}, while the three
 dimensional prognostic variables are solved with a longer time step that is a
 multiple of \np{rn\_rdt}. The trend in the barotropic momentum due to bottom
+equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three
+dimensional prognostic variables are solved with the longer time step
+of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom
 friction appropriate to this method is that given by the selected parameterisation
 ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities
 …
 \end{enumerate}
 Note that the use of an implicit formulation
+Note that the use of an implicit formulation within the barotropic loop
 for the bottom friction trend means that any limiting of the bottom friction coefficient
 in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time
 splitting. This is because the major contribution to bottom friction is likely to come from
+the barotropic component which uses the unrestricted value of the coefficient.
+The implicit formulation takes the form:
+the barotropic component which uses the unrestricted value of the coefficient. However, if the
+limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas
+applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} )
+which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}.
+Otherwise, the implicit formulation takes the form:
 \begin{equation} \label{Eq_zdfbfr_implicitts}
  \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ]
 …
 The essential goal of the parameterization is to represent the momentum
 exchange between the barotropic tides and the unrepresented internal waves
 induced by the tidal ßow over rough topography in a stratified ocean.
+induced by the tidal flow over rough topography in a stratified ocean.
 In the current version of \NEMO, the map is built from the output of
 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}.

trunk/DOC/TexFiles/Chapters/Introduction.tex

Property svn:executable deleted

-                      r2570
+                      r3294
 \citep{OASIS2006}. Two-way nesting is also available through an interface to the
 AGRIF package (Adaptative Grid Refinement in \textsc{Fortran}) \citep{Debreu_al_CG2008}.
+The interface code for coupling to an alternative sea ice model (CICE, \citet{Hunke2008}) is now
+available although this is currently only designed for global domains, without the use of AGRIF.
 Other model characteristics are the lateral boundary conditions (chapter~\ref{LBC}).
 …
 concern of improving the model performance.
+ \vspace{1cm}
+$\bullet$ The main modifications from NEMO/OPA v3.3 and  v3.4 are :\\
+\begin{enumerate}
+\item finalisation of above iso-neutral mixing \citep{Griffies_al_JPO98}";
+\item "Neptune effect" parametrisation;
+\item horizontal pressure gradient suitable for s-coordinate;
+\item semi-implicit bottom friction;
+\item finalisation of the merge of passive and active tracers advection-diffusion modules;
+\item a new bulk formulae (so-called MFS);
+\item use fldread for the off-line tracer component (OFF\_SRC) ;
+\item use MPI point to point communications  for north fold;
+\item diagnostic of transport ;
+\end{enumerate}

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

Context Navigation

Legend:

trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex

trunk/DOC/TexFiles/Chapters/Annex_A.tex

trunk/DOC/TexFiles/Chapters/Annex_B.tex

trunk/DOC/TexFiles/Chapters/Annex_C.tex

trunk/DOC/TexFiles/Chapters/Annex_D.tex

trunk/DOC/TexFiles/Chapters/Annex_E.tex

trunk/DOC/TexFiles/Chapters/Annex_ISO.tex

trunk/DOC/TexFiles/Chapters/Chap_ASM.tex

trunk/DOC/TexFiles/Chapters/Chap_CFG.tex

trunk/DOC/TexFiles/Chapters/Chap_Conservation.tex

trunk/DOC/TexFiles/Chapters/Chap_DIA.tex

trunk/DOC/TexFiles/Chapters/Chap_DOM.tex

trunk/DOC/TexFiles/Chapters/Chap_DYN.tex

trunk/DOC/TexFiles/Chapters/Chap_LBC.tex

trunk/DOC/TexFiles/Chapters/Chap_LDF.tex

trunk/DOC/TexFiles/Chapters/Chap_MISC.tex

trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex

trunk/DOC/TexFiles/Chapters/Chap_Model_Basics_zstar.tex

trunk/DOC/TexFiles/Chapters/Chap_OBS.tex

trunk/DOC/TexFiles/Chapters/Chap_SBC.tex

trunk/DOC/TexFiles/Chapters/Chap_STP.tex

trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex

trunk/DOC/TexFiles/Chapters/Introduction.tex

Download in other formats: