Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex	(revision 3236)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex	(revision 3267)
@@ -2,12 +2,67 @@
 % Iso-neutral diffusion :
 % ================================================================
-\chapter{Griffies's iso-neutral diffusion}
+\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}
+\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 define a scheme inspired by \citet{Griffies_al_JPO98}, but formulated within the \NEMO
+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.
 
@@ -613,11 +668,11 @@
 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 \nlv{ln\_botmix\_grif=.false.}, but left unmasked,
-giving bottom mixing, if \nlv{ln\_botmix\_grif=.true.}.
-
-The default option \nlv{ln\_botmix\_grif=.false.} is suitable when the
-bbl mixing option is enabled (\key{trabbl}, with \nlv{nn\_bbl\_ldf=1}),
+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, \nlv{ln\_botmix\_grif=.true.} may be necessary.
+bbl mixing, \np{ln\_botmix\_grif}=true may be necessary.
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[h] \begin{center}
@@ -636,6 +691,6 @@
       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 \smnlv{ln\_botmix\_grif=.false.}, but left
-      unmasked, giving bottom mixing, if \smnlv{ln\_botmix\_grif=.true.}}
+      line) are masked if \np{botmix\_grif}=.false., but left
+      unmasked, giving bottom mixing, if \np{botmix\_grif}=.true.}
  \end{center} \end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -665,5 +720,6 @@
 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}
+\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
@@ -671,5 +727,5 @@
 \subsubsection{Linear slope tapering within the surface mixed layer}\label{sec:triad:lintaper}
 This is the option activated by the default choice
-\nlv{ln\_triad\_iso=.false.}. Slopes $\tilde{r}_i$ relative to
+\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
@@ -798,7 +854,9 @@
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-\subsubsection{Additional truncation of skew iso-neutral flux components}
-The alternative option is activated by setting \nlv{ln\_triad\_iso =
-  .true.}. This retains the same tapered slope $\rML$  described above for the
+\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
@@ -839,9 +897,9 @@
 % Skew flux formulation for Eddy Induced Velocity :
 % ================================================================
-\section{Eddy induced advection and its formulation as a skew flux}
+\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 (\key{traldf\_eiv} defined),
+ 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-
@@ -852,22 +910,36 @@
 
 The eddy induced velocity is given by:
-\begin{equation} \label{eq:triad:eiv_v}
+\begin{subequations} \label{eq:triad:eiv}
+\begin{equation}\label{eq:triad:eiv_v}
 \begin{split}
- u^* & = - \frac{1}{e_{3}}\;          \partial_k \left( A_{e} \; \tilde{r}_1 \right)   \\
- v^* & = - \frac{1}{e_{3}}\;          \partial_k \left( A_{e} \; \tilde{r}_2 \right)   \\
-w^* & =    \frac{1}{e_{1}e_{2}}\; \left\{ \partial_i  \left( e_{2} \, A_{e} \; \tilde{r}_1 \right)
-							    + \partial_j  \left( e_{1} \, A_{e} \;\tilde{r}_2 \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 $\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 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
+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
@@ -883,16 +955,16 @@
 &=
 \begin{pmatrix}
- 	        { - \partial_k \left( e_{2} \, A_{e} \; \tilde{r}_1 \right) \; T \;} 	 	\\
- 		{+ \partial_i  \left( e_{2} \, A_{e} \; \tilde{r}_1 \right) \; T \;}	 \\
+ 	        { - \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} \; \tilde{r}_1  \; T \right) \;}  \\
- 		{+ \partial_i  \left( e_{2} \, A_{e} \; \tilde{r}_1  \; T \right) \;}	 \\
+ 	        { - \partial_k \left( e_{2} \, \psi_1  \; T \right) \;}  \\
+ 		{+ \partial_i  \left( e_{2} \,\psi_1 \; T \right) \;}	 \\
 \end{pmatrix}
  +
 \begin{pmatrix}
- 	        {+ e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_k T}  \\
- 		{ - e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_i  T}	 \\
+ 	        {+ e_{2} \, \psi_1  \; \partial_k T}  \\
+ 		{ - e_{2} \, \psi_1  \; \partial_i  T}	 \\
 \end{pmatrix}
 \end{split}
@@ -902,6 +974,6 @@
 \begin{equation} \label{eq:triad:eiv_skew_ijk}
 \textbf{F}_\mathrm{eiv}^T = \begin{pmatrix}
- 	        {+ e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_k T}   \\
- 		{ - e_{2} \, A_{e} \; \tilde{r}_1  \; \partial_i  T}	 \\
+ 	        {+ e_{2} \, \psi_1  \; \partial_k T}   \\
+ 		{ - e_{2} \, \psi_1  \; \partial_i  T}	 \\
                                  \end{pmatrix}
 \end{equation}
@@ -909,22 +981,22 @@
 \begin{equation}\label{eq:triad:eiv_skew_physical}
 \begin{split}
- f^*_1 & = \frac{1}{e_{3}}\; A_{e} \; \tilde{r}_1 \partial_k T   \\
- f^*_2 & = \frac{1}{e_{3}}\; A_{e} \; \tilde{r}_2 \partial_k T   \\
- f^*_3 & =  -\frac{1}{e_{1}e_{2}}\; A_{e} \left\{ e_{2} \tilde{r}_1 \partial_i T
-   + e_{1} \tilde{r}_2 \partial_j T \right\}. \\
+ 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$ are relative to geopotentials.
+$\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  A_{e} \; \tilde{r}_1 \partial_k T\right)
-  + \frac{\partial}{\partial j} \left( e_1  A_{e} \;
-    \tilde{r}_2 \partial_k T\right)
- -  \frac{\partial}{\partial k} A_{e} \left( e_{2} \tilde{r}_1 \partial_i T
-   + e_{1} \tilde{r}_2 \partial_j T \right)  \right]
+  \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
@@ -976,5 +1048,7 @@
 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
+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
@@ -1064,5 +1138,5 @@
 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 \nlv{ln\_botmix\_grif} has
+$u$-point is masked. The namelist parameter \np{ln\_botmix\_grif} has
 no effect on the eddy-induced skew-fluxes.
 
@@ -1079,5 +1153,5 @@
 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 \nlv{ln\_triad\_iso}.
+unaffected by the value of \np{ln\_triad\_iso}.
 
 The justification for this linear slope tapering is that, for $A_e$
@@ -1094,5 +1168,5 @@
 
 \subsection{Streamfunction diagnostics}\label{sec:triad:sfdiag}
-Where the namelist parameter \nlv{ln\_botmix\_grif=.true.}, diagnosed
+Where the namelist parameter \np{ln\_botmix\_grif}=true, diagnosed
 mean eddy-induced velocities are output. Each time step,
 streamfunctions are calculated in the $i$-$k$ and $j$-$k$ planes at
@@ -1104,83 +1178,18 @@
 \begin{equation}
   \label{eq:triad:sfdiagi}
-  {\psi_{[i]}}_{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}
-
-\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}_\mathrm{eiv}(T) \; T \;dv  \equiv 0
-\end{align*}
-The discrete form of its left hand side is obtained using \eqref{eq:triad:allskewflux}
-\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-luatex.tex"
-%%% End:
+  {\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}
+\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}
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3236)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3267)
@@ -123,10 +123,10 @@
 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)		\\
@@ -139,7 +139,7 @@
 \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}
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_LDF.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 3236)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_LDF.tex	(revision 3267)
@@ -21,5 +21,8 @@
 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--------------------------------------------
@@ -128,10 +131,10 @@
 $\ $\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 
@@ -275,5 +278,6 @@
 
 \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 
@@ -332,20 +336,24 @@
 \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, 
Index: branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 3236)
+++ branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 3267)
@@ -491,5 +491,6 @@
 \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:
@@ -536,7 +537,7 @@
 (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 
