Changeset 3218

Timestamp:

2011-12-13T11:31:24+01:00 (12 years ago)

Author:

agn

Message:

agn: changes to documentation in dev_NEMO_MERGE_2011

--Reordered NEMO_book.tex, changed from 11pt->12pt,

only needs pstricks for logos

--more abbreviations in mathabbrev.sty
--Corrected bug in Abstracts_Foreword.tex that put whole

text into small

--Clarified discussion of LDF operator rotation in
Chap_Model_Basics.tex, Chap_LDF.tex, Annex_B.tex
--Discuss surface tapering of Griffies triads, how they

work in s-coordinates, eddy-induced velocities in
Annex_ISO.tex. Extra ref in Biblio.bib
Extra figures Fig_bdry_triads.pdf, Fig_triad_MLB.pdf
Modified namelist section namtra_ldf

--Minor corrections to Chap_DIA.tex, Chap_CFG.tex,

Chap_SBC.tex, Annex_A.tex, Chap_DOM.tex, Chap_ASM.tex

Location:

branches/2011/dev_NEMO_MERGE_2011/DOC

Files:

• NEMO_book.tex (modified) (15 diffs)
• TexFiles/Biblio/Biblio.bib (modified) (4 diffs)
• TexFiles/Chapters/Abstracts_Foreword.tex (modified) (1 diff)
• TexFiles/Chapters/Annex_A.tex (modified) (7 diffs)
• TexFiles/Chapters/Annex_B.tex (modified) (11 diffs)
• TexFiles/Chapters/Annex_ISO.tex (modified) (27 diffs)
• TexFiles/Chapters/Chap_ASM.tex (modified) (1 diff)
• TexFiles/Chapters/Chap_CFG.tex (modified) (2 diffs)
• TexFiles/Chapters/Chap_DIA.tex (modified) (2 diffs)
• TexFiles/Chapters/Chap_DOM.tex (modified) (3 diffs)
• TexFiles/Chapters/Chap_LDF.tex (modified) (2 diffs)
• TexFiles/Chapters/Chap_Model_Basics.tex (modified) (7 diffs)
• TexFiles/Chapters/Chap_SBC.tex (modified) (4 diffs)
• TexFiles/Figures/Fig_bdry_triads.pdf (added)
• TexFiles/Figures/Fig_triad_MLB.pdf (added)
• TexFiles/Namelist/namtra_ldf (modified) (1 diff)
• TexFiles/math_abbrev.sty (modified) (1 diff)

branches/2011/dev_NEMO_MERGE_2011/DOC/NEMO_book.tex

@@ -4 +4 @@
 % (C) Xavier Perseguers 2002 - xavier.perseguers@epfl.ch
 
-\documentclass[a4paper,12pt]{book}
+\documentclass[a4paper,11pt]{book}
 
 %  makeindex NEMO_book     <== to regenerate the index
@@ -14 +14 @@
 
 \usepackage[french]{babel}
-\usepackage{fancyhdr}
-
-\usepackage{alltt}      %%  alltt for namelist
-\usepackage{verbatim}   %%  alltt for namelist
-
-%hyperref 
+%\usepackage{color}
+\usepackage{xcolor}
+%\usepackage{graphics}           % allows insertion of pictures
+\usepackage{graphicx}            % allows insertion of
+                                        % pictures
+\usepackage[capbesideposition={top,center}]{floatrow}
+\floatsetup[table]{style=plaintop}
+\usepackage[margin=10pt,font={small},labelsep=colon,labelfont={bf}]{caption}
+%\usepackage{colortbl}                           % gives coloured panels behind table columns
+
+%hyperref
 \usepackage[               %
   pdftitle={NEMO ocean engine},  %
@@ -31 +36 @@
   linkcolor=blue,anchorcolor=blue,  %
   citecolor=blue,filecolor=blue,    %
-  menucolor=blue,pagecolor=blue, %
+ menucolor=blue,                    %
   urlcolor=blue]{hyperref}
 %  usage of exteranl hyperlink :  \href{mailto:my_address@wikibooks.org}{my\_address@wikibooks.org}
@@ -39 +44 @@
 
 
-
-%\usepackage{amssymb}
+%%%% page styles etc................
+\usepackage{fancyhdr}
 \pagestyle{fancy}
 % with this we ensure that the chapter and section
@@ -59 +64 @@
 }
 
+
+%%%%  Section number in Margin.......
+% typeset the number of each section in the left margin, with the start of each instance of
+% sectional heading text aligned with the left hand edge of  the body text.
+\makeatletter
+\def\@seccntformat#1{\protect\makebox[0pt][r]{\csname the#1\endcsname\quad}}
+\makeatother
+
 % Leave blank pages completely empty, w/o header
 \makeatletter
@@ -70 +83 @@
 \makeatother
 
-\usepackage{makeidx}
-
-\usepackage{minitoc}          %in french : \usepackage[french]{minitoc}
-
+%%%% define the chapter  style ................
+\usepackage{minitoc}          %In French : \usepackage[french]{minitoc}
 %\usepackage{mtcoff}          % invalidate the use of minitocs
-\usepackage[latin1]{inputenc}
-\usepackage{graphics}            % allows insertion of pictures
-\usepackage{times}               % use vector fonts instead of bitmap
 \usepackage{fancybox}
-\usepackage{graphicx}
-%\usepackage{color}
-%\usepackage{colortbl}
-
-%%%% add by Word2Tex
-%\documentclass{amsart}
-\usepackage{latexsym}
-\usepackage{amssymb}
-\usepackage{amsmath}
-\usepackage{graphicx}
-\allowdisplaybreaks[1]           % allow page breaks in the middle of equations
-%%%%
-
-%%%%
-%\usepackage{mathtime}      % font for illustrator to work (belleek fonts )
-%%%%
-
-%\usepackage{showidx}            % show the index entry
-
-
-
-%%%% define the style of title for chapter  %%%
-\usepackage{pstricks}
+
 \makeatletter
 \def\LigneVerticale{\vrule height 5cm depth 2cm\hspace{0.1cm}\relax}
@@ -108 +94 @@
 \def\GrosCarreAvecUnChiffre#1{%
   \rlap{\vrule height 0.8cm width 1cm depth 0.2cm}%
-  \rlap{\hbox to 1cm{\hss\mbox{\white #1}\hss}}%
+ \rlap{\hbox to 1cm{\hss\mbox{\color{white} #1}\hss}}%
   \vrule height 0pt width 1cm depth 0pt}
 \def\GrosCarreAvecTroisChiffre#1{%
   \rlap{\vrule height 0.8cm width 1.6cm depth 0.2cm}%
-  \rlap{\hbox to 1.5cm{\hss\mbox{\white #1}\hss}}%
+ \rlap{\hbox to 1.5cm{\hss\mbox{\color{white} #1}\hss}}%
   \vrule height 0pt width 1cm depth 0pt}
 
 \def\@makechapterhead#1{\hbox{%
-    \huge 
+   \huge
     \LignesVerticales
     \hspace{-0.5cm}%
@@ -126 +112 @@
 }\par\vskip 1cm}
 \def\@makeschapterhead#1{\hbox{%
-    \huge 
+   \huge
     \LignesVerticales
     %\hspace{0.5cm}%
     \hbox{#1}%
 }\par\vskip 2cm}
+\makeatother
 
 %\def\thechapter{\Roman{chapter}}      % chapter number to be Roman
 
-%%%% end style chapter %%%
-
-%%% use maths shortcuts :
-%\newcommand{\vect} [1] { \rm{\textbf{#1}} }    % vector style: non-italic bold
-\usepackage{./TexFiles/math_abbrev}
+
+%%%%           Mathematics...............
+%\documentclass{amsart}
+\usepackage{xspace}                              % helpd ensure correct spacing after macros
+\usepackage{latexsym}
+\usepackage{amssymb}
+\usepackage{amsmath}
+\allowdisplaybreaks[1]           % allow page breaks in the middle of equations
+\usepackage{./TexFiles/math_abbrev}    % use maths shortcuts
+
+
+\usepackage{times}                % use times font for text
+%\usepackage{mathtime}                          % font for illustrator to work (belleek fonts )
+%\usepackage[latin1]{inputenc}                % allows some unicode removed (agn)
 
 
 %%% essai commande
-\newcommand{\nl} [1] {\texttt{\small {\textcolor{blue}{#1}} } } 
-
-%%% index commands :   
+\newcommand{\nl} [1] {\texttt{\small {\textcolor{blue}{#1}} } }
+\newcommand{\nlv} [1] {\texttt{\footnotesize#1}\xspace}
+\newcommand{\smnlv} [1] {\texttt{\scriptsize#1}\xspace}
+
+%%%% namelist & code display................................
+\usepackage{alltt}      %%  alltt for namelist
+\usepackage{verbatim}   %%  alltt for namelist
+% namelists
+\newcommand{\namdisplay} [1] {
+\begin{alltt}
+{\tiny \verbatiminput{./TexFiles/Namelist/#1}}
+\end{alltt}
+  \vspace{-10pt}
+}
+% code display
+%\newcommand{\codedisplay} [1] { \begin{alltt} {\tiny  {\begin{verbatim} {#1}} \end{verbatim} }  \end{alltt}   }
+
+
+
+%%%% commands for working with text................................
+% command to "comment out" portions of text ({} argument) or not ({#1} argument)
+\newcommand{\amtcomment}[1]{}    % command to "commented out" portions of text or not (#1 in argument)
+\newcommand{\sgacomment}[1]{}    % command to "commented out" portions of
+\newcommand{\gmcomment}[1]{}     % command to "commented out" portions of
+%                                               % text that span line breaks
+%Red (NR) or Yellow(WARN)
+%\newcommand{\NR} {\colorbox{red}{#1}}
+%\newcommand{\WARN} {{ \colorbox{yellow}{#1}} }
+
+
+
+%%% index commands......................
+\usepackage{makeidx}
+%\usepackage{showidx}            % show the index entry
+
 \newcommand{\mdl} [1] {\textit{#1.F90}\index{Modules!#1}}         %module (mdl)
 \newcommand{\rou} [1] {\textit{#1}\index{Routines!#1}}            %module (routine)
@@ -153 +181 @@
 \newcommand{\ifile} [1] {\textit{#1.nc}\index{Input NetCDF files!#1.nc}}   %input NetCDF files (.nc)
 \newcommand{\key} [1] {\textbf{key\_#1}\index{CPP keys!key\_#1}}  %key_cpp (key)
-\newcommand{\NEMO} {\textit{NEMO }}                %NEMO (nemo)
-% command to "commented out" portions of text ({} argument) or not ({#1} argument)
-\newcommand{\amtcomment}[1]{}    % command to "commented out" portions of text or not (#1 in argument)
-\newcommand{\sgacomment}[1]{}    % command to "commented out" portions of 
-\newcommand{\gmcomment}[1]{}     % command to "commented out" portions of 
-%                                               % text that span line breaks
-\newcommand{\alpbet} {\left(\alpha / \beta \right)}   % alpha/beta  for slp computation
-
-%--------------------------------------------namlist---------------------------------------------------------
-\newcommand{\namdisplay} [1] { 
-   \begin{alltt}  {{\tiny \verbatiminput{./TexFiles/Namelist/#1}}}  \end{alltt} 
-                               \vspace{-10pt}   }
-%--------------------------------------------code display---------------------------------------------------------
-%\newcommand{\codedisplay} [1] { \begin{alltt} {\tiny  {\begin{verbatim} {#1}} \end{verbatim} }  \end{alltt}   }
-%--------------------------------------------------------------------------------------------------------------
-%Red (NR) or Yellow(WARN)
-%\newcommand{\NR} {\colorbox{red}{#1}}
-%\newcommand{\WARN} {{ \colorbox{yellow}{#1}} }
-
-%---------------------%%%%%  Section number in Margin  %%%%%-----------------------------
-% typeset the number of each section in the left margin, with the start of each instance of
-% sectional heading text aligned with the left hand edge of  the body text.
-\makeatletter 
-\def\@seccntformat#1{\protect\makebox[0pt][r]{\csname the#1\endcsname\quad}} 
-\makeatother 
-%--------------------------------------------------------------------------------------------------------------
-
-%----------------------------%%%%%   Bilbiography   %%%%%%-------------------------------%
+\newcommand{\NEMO} {\textit{NEMO}\xspace}                %NEMO (nemo)
+
+%%%%   Bibliography   .............
 \usepackage[nottoc, notlof, notlot]{tocbibind}
 \usepackage[square, comma]{natbib}
 \bibpunct{[}{]}{,}{a}{}{;}                           %suppress "," after "et al."
 \providecommand{\bibfont}{\small}
-%--------------------------------------------------------------------------------------------------------------
+
 
 % ================================================================
@@ -191 +194 @@
 % ================================================================
 
-\title{ 
+\usepackage{pstricks}
+\title{
 \psset{unit=1.1in,linewidth=4pt}    %parameters of the units for pstricks
 \rput(-1.4,2){ \includegraphics[width=0.2\textwidth]{./TexFiles/Figures/logo_CNRS.pdf}           }
@@ -205 +209 @@
 \rule{345pt}{1.5pt} \\
 \vspace{0.45cm}
- {\Huge NEMO ocean engine} 
-\rule{345pt}{1.5pt} \\ 
+{\Huge NEMO ocean engine}
+\rule{345pt}{1.5pt} \\
  }
 %{ -- Draft --}   }
@@ -220 +224 @@
 
 
-%%\includegraphics[width=0.2\textwidth]{./TexFiles/Figures/logo_NERC.pdf} \\
-
-%%    \hbox{\mbox{%
-%%       \hspace{4pt}%
-%%         \fbox{\includegraphics[width=3em]{./TexFiles/Figures/logo_CNRS.pdf}}%
-%%        \hspace{4pt}
-%%        }  }%
-
-\author{ 
+\author{
 \Large Gurvan Madec, and the NEMO team  \\
  \texttt{\small gurvan.madec@locean-ipsl.umpc.fr} \\
@@ -257 +253 @@
 \frontmatter
 
-\tableofcontents              % generate a table of content
-%\listoffigures               % generate a table of content
-%\listoftables             % generate a table of content
+\tableofcontents              % generate a table of contents
+%\listoffigures               % generate a list  of figures
+%\listoftables                     % generate a list of tables
 
 \mainmatter
@@ -283 +279 @@
 \include{./TexFiles/Chapters/Chap_STP}       % Time discretisation (time stepping strategy)
 
-\include{./TexFiles/Chapters/Chap_DOM}       % Space discretisation 
+\include{./TexFiles/Chapters/Chap_DOM}       % Space discretisation
 
 \include{./TexFiles/Chapters/Chap_TRA}       % Tracer advection/diffusion equation
@@ -318 +314 @@
 \include{./TexFiles/Chapters/Annex_C}        % Discrete invariants of the eqs.
 \include{./TexFiles/Chapters/Annex_D}        % Coding rules
-\include{./TexFiles/Chapters/Annex_ISO}                     % Isoneutral diffusion using triads 
-%\include{./TexFiles/Chapters/Annex_E}                   % Notes on some on going staff (no included in the DOC) 
-%\include{./TexFiles/Chapters/Annex_Fox-Kemper}   % Notes on Fox-Kemper (no included in the DOC) 
-%\include{./TexFiles/Chapters/Annex_EVP}           % Notes on EVP (no included in the DOC) 
+\include{./TexFiles/Chapters/Annex_ISO}                     % Isoneutral diffusion using triads
+%\include{./TexFiles/Chapters/Annex_E}                   % Notes on some on going staff (no included in the DOC)
+%\include{./TexFiles/Chapters/Annex_Fox-Kemper}   % Notes on Fox-Kemper (no included in the DOC)
+%\include{./TexFiles/Chapters/Annex_EVP}           % Notes on EVP (no included in the DOC)
 
 % ================================================================

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Biblio/Biblio.bib

@@ -705 +705 @@
 }
 
+
+@ARTICLE{Danabasoglu_al_2008,
+   author = {G. Danabasoglu and R. Ferrari and J. C. McWilliams},
+   title = {Sensitivity of an ocean general circulation model to a parameterization of near-surface eddy fluxes},
+   volume = {21},
+   year = {2008},
+   doi = {10.1175/2007JCLI1508.1},
+   journal = {J. Climate},
+   pages = {1192--1208},
+        url = {http://dx.doi.org/10.1175/2007JCLI1508.1}
+}
 @ARTICLE{Dandonneau_al_S04,
   author = {Y. Dandonneau and C. Menkes and T. Gorgues and G. Madec},
@@ -1059 +1070 @@
   pages = {14703--14726}
 }
-
+@ARTICLE{Gerdes1991,
+   Author = {Gerdes, R{\"u}diger and K{\"o}berle, Cornelia and Willebrand, J{\"u}rgen},
+   Doi = {10.1007/BF00210006},
+   Journal = {Clim. Dynamics},
+   Number = {4},
+   Pages = {211--226},
+   Title = {The influence of numerical advection schemes on the results of ocean general circulation models},
+   Url = {http://dx.doi.org/10.1007/BF00210006},
+   Volume = {5},
+   Year = {1991},
+}
 @TECHREPORT{Gibson_TR86,
   author = {J. K. Gibson},
@@ -1133 +1154 @@
   pages = {1--46},
   doi = {10.1016/j.ocemod.2008.08.007},
-  url = {http://dx.doi.org/}
+  url = {http://dx.doi.org/10.1016/j.ocemod.2008.08.007}
 }
 
@@ -2301 +2322 @@
   journal = JAS,
   year = {1972},
-  volume = {29}
+  volume = {29},
   pages = {959--970}
 }

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Abstracts_Foreword.tex

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

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_A.tex

@@ -13 +13 @@
 % 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 
@@ -62 +62 @@
 % continuity equation
 % ================================================================
-\section{Continuity Equation in $s-$coordinate}
+\section{Continuity Equation in $s-$coordinates}
 \label{Apdx_A_continuity}
 
@@ -128 +128 @@
 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    \\
@@ -429 +429 @@
 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 
@@ -451 +452 @@
 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}
@@ -495 +496 @@
 \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
@@ -502 +503 @@
 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. 
 
 

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_B.tex

@@ -16 +16 @@
 \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}
  1 \hfill & 0 \hfill & {-\sigma _1 } \hfill \\
  0 \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)
@@ -43 +44 @@
 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 \\
@@ -120 +122 @@
 \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 \\
@@ -129 +131 @@
 \end{align*}
 \end{subequations}
- 
+\addtocounter{equation}{-2}
 
 % ================================================================
@@ -137 +139 @@
 \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}
@@ -149 +152 @@
  {-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}
  1 \hfill & 0 \hfill & {-a_1 } \hfill \\
  0 \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}
+ 1 \hfill & 0 \hfill & {-r _1 } \hfill \\
+ 0 \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}
+ 1 \hfill & 0 \hfill &0 \hfill \\
+ 0 \hfill & 1 \hfill & 0 \hfill \\
+0 \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.
+
 
 % ================================================================
@@ -211 +289 @@
 \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)    \\
@@ -224 +302 @@
  {\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)
@@ -249 +327 @@
 \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*}
@@ -255 +333 @@
 \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}
@@ -280 +358 @@
 \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.

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Annex_ISO.tex

@@ -1 +1 @@
 % ================================================================
-% Iso-neutral diffusion : 
+% Iso-neutral diffusion :
 % ================================================================
 \chapter{Griffies's iso-neutral diffusion}
-\label{Apdx_C}
+\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
+\label{sec:triad:iso}
+
+We define 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}
+      1&0&-r_1\mystrut \\
+      0&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}
@@ -33 +57 @@
     }{\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*}
@@ -76 +123 @@
 \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
@@ -89 +137 @@
 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.
@@ -95 +143 @@
 \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 ] }.
@@ -160 +207 @@
 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
@@ -202 +248 @@
 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} }
@@ -227 +273 @@
 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}
@@ -234 +280 @@
 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} }
@@ -260 +306 @@
   \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)$
@@ -270 +316 @@
 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}
@@ -284 +330 @@
   \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+
@@ -311 +357 @@
   &= -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(
@@ -343 +390 @@
 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}
@@ -350 +397 @@
 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} }
@@ -358 +405 @@
   \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
@@ -365 +412 @@
 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
@@ -376 +423 @@
 \[
 \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)
 \]
@@ -390 +437 @@
 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[
@@ -427 +513 @@
 \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}
@@ -437 +523 @@
 \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}
@@ -454 +541 @@
 \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}
@@ -466 +553 @@
 \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}
@@ -486 +573 @@
   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} }
@@ -506 +593 @@
 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 \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}),
+or  for simple idealized  problems. For setups with topography without
+bbl mixing, \nlv{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 \smnlv{ln\_botmix\_grif=.false.}, but left
+      unmasked, giving bottom mixing, if \smnlv{ln\_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}
+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
+\nlv{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}
+ 1 \hfill & 0 \hfill & {-\rML[1]}\hfill \\
+ 0 \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}
+The alternative option is activated by setting \nlv{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}
+ 1 \hfill & 0 \hfill & {-\rML[1]^*}\hfill \\
+ 0 \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 and its formulation as a 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),
+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{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_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\}
 \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 $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
+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} \, A_{e} \; \tilde{r}_1 \right) \; T \;}     \\
+      {+ \partial_i  \left( e_{2} \, A_{e} \; \tilde{r}_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} \, A_{e} \; \tilde{r}_1  \; T \right) \;}  \\
+      {+ \partial_i  \left( e_{2} \, A_{e} \; \tilde{r}_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}   \\
+\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} \, A_{e} \; \tilde{r}_1  \; \partial_k T}   \\
+      { - e_{2} \, A_{e} \; \tilde{r}_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.
-
-
-
+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}}\; 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\}. \\
+\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.
+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]
+\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 \nlv{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 \nlv{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 \nlv{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
+$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_{[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
@@ -608 +1116 @@
 
 
-Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane. 
+Demonstration for the conservation of the tracer variance in the (\textbf{i},\textbf{j}) plane.
 
 This have to be moved in an Appendix.
@@ -614 +1122 @@
 The continuous property to be demonstrated is :
 \begin{align*}
-\int_D \nabla \cdot \textbf{F}_{eiv}(T) \; T \;dv  \equiv 0
+\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_eiv_skew}
+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}]         
+ \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\}   
+- \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}]         
+  &\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\}       
+- &\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} 
+ \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\ \ \ \:}      
+&&+{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} 
+&&+{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\ \ \ \:}       
+&&+{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\ \ \ \:} 
+&&-{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}] 
+&&-{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{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 
+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: 
+%%% Local Variables:
+%%% TeX-master: "../../NEMO_book-luatex.tex"
+%%% End:

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_ASM.tex

@@ -100 +100 @@
  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{aligned} \right.,
 \end{equation}
 where

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_CFG.tex

@@ -219 +219 @@
 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}. 
@@ -227 +227 @@
 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.

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DIA.tex

@@ -236 +236 @@
 \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.
 
@@ -321 +321 @@
 \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 

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_DOM.tex

@@ -81 +81 @@
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{table}[!tb] \label{Tab_cell}
+\begin{table}[!tb]
 \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
 \hline
@@ -93 +93 @@
 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. 
@@ -851 +851 @@
 \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}

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_LDF.tex

@@ -220 +220 @@
 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 
@@ -369 +369 @@
 
 
-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 

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_Model_Basics.tex

@@ -1094 +1094 @@
 % 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}
 
@@ -1134 +1134 @@
 
 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 
@@ -1157 +1157 @@
 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}):
@@ -1186 +1186 @@
 
 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}
@@ -1213 +1217 @@
 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}
@@ -1227 +1231 @@
 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:
@@ -1239 +1243 @@
 
 
-\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 

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Chapters/Chap_SBC.tex

@@ -246 +246 @@
 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}
@@ -619 +620 @@
 \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})
+\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
 \end{itemize}
 % -------------------------------------------------------------------------------------------------------------
@@ -862 +863 @@
 %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:
 
-}
+%}
 
 % ================================================================
@@ -1078 +1079 @@
 \begin{description}
 
-In order to read a neutral drag coeff, from an external data source (i.e. a wave model), the 
+\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.}$. 

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/Namelist/namtra_ldf

@@ -1 +1 @@
 !-----------------------------------------------------------------------
-&namtra_ldf    !   lateral diffusion scheme for tracer 
+&namtra_ldf    !   lateral diffusion scheme for tracer
 !-----------------------------------------------------------------------
-   !                       !  Type of the operator : 
-   ln_traldf_lap    =  .true.   !  laplacian operator       
-   ln_traldf_bilap  =  .false.  !  bilaplacian operator     
+   !                       !  Type of the operator :
+   ln_traldf_lap    =  .true.   !  laplacian operator
+   ln_traldf_bilap  =  .false.  !  bilaplacian operator
    !                       !  Direction of action  :
-   ln_traldf_level  =  .false.  !  iso-level                
-   ln_traldf_hor    =  .false.  !  horizontal (geopotential)            (require "key_ldfslp" when ln_sco=T)
-   ln_traldf_iso    =  .true.   !  iso-neutral                          (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !  griffies skew flux formulation       (require "key_ldfslp")
-   ln_traldf_gdia   =  .false.  !  griffies operator strfn diagnostics  (require "key_ldfslp")
-   ln_triad_iso     =  .false.  !  griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp")
-   ln_botmix_grif   =  .false.  !  griffies operator with lateral mixing on bottom (require "key_ldfslp")
-   !                       !  Coefficient
+   ln_traldf_level  =  .false.  !  iso-level
+   ln_traldf_hor    =  .false.  !  horizontal (geopotential)            (requires "key_ldfslp" when ln_sco=T)
+   ln_traldf_iso    =  .true.   !  iso-neutral                          (requires "key_ldfslp")
+   !                 !  Griffies parameters
+   ln_traldf_grif   =  .false.  !  use griffies triad formulation       (requires "key_ldfslp")
+   ln_traldf_gdia   =  .false.  !  output griffies strfn diagnostics    (requires "key_ldfslp")
+   ln_triad_iso     =  .false.  !  isoneutral diff'n triads => pure lateral mixing in ML (requires "key_ldfslp")
+   ln_botmix_grif   =  .false.  !  griffies operator with lateral mixing on bottom (requires "key_ldfslp")
+   !                       !  Coefficients
+   ! rn_aeiv_0 is ignored with Held-Larichev spatially varying aeiv (key_traldf_c2d & key_orca_r2, _r1 or _r05)
+   rn_aeiv_0        =  2000.    !  eddy induced velocity coefficient [m2/s]    (requires "key_traldf_eiv")
    rn_aht_0         =  2000.    !  horizontal eddy diffusivity for tracers [m2/s]
-   rn_ahtb_0        =     0.    !  background eddy diffusivity for ldf_iso [m2/s]
-   rn_aeiv_0        =  2000.    !  eddy induced velocity coefficient [m2/s]    (require "key_traldf_eiv")
+   rn_ahtb_0        =     0.    !  background eddy diffusivity for ldf_iso [m2/s] (normally=0; not used with Griffies)
 /

branches/2011/dev_NEMO_MERGE_2011/DOC/TexFiles/math_abbrev.sty

@@ -6 +6 @@
 \def\stwelfth{\sfrac{1}{12}}
 \def\sthirtysixth{\sfrac{1}{36}}
-\def\grad{\boldsymbol{\nabla}}
-\def\gradh{\boldsymbol{\nabla}_h}
+% \newcommand{\grad}{\boldsymbol{\nabla}}
+% \newcommand{\gradh}{\boldsymbol{\nabla}_h}
+\newcommand{\grad}{\nabla}
+\newcommand{\gradh}{\nabla_h}
 \def\bgamma\boldsymbol{\gamma}
 \def\rdt{\Delta t}
 \newcommand{\ew}[3]{{e_{3#1}}_{\,#2}^{\,#3} }
-\newcommand{\vect}[1]{ \rm{\textbf{#1}} }    % vector style: non-italic bold
-\def\deg{\degres}                            % degrees  (NB: \r{} can also be used)
+\newcommand{\vect}[1]{\ensuremath{\mathbf{#1}}} % for vectors                                % non-italic bold
+\newcommand{\Div}{\grad\cdot}
+% \newcommand{\curl}{\boldsymbol{\nabla} \times} % for curl
+\newcommand{\curl}{\nabla \times} % for curl
+\newcommand{\pd}[2][]{\frac{\partial #1}{\partial #2}}
+\def\deg{\degres}                            % degrees  (NB: \r{} can % % also be used)
+\newcommand{\alpbet} {\left(\alpha / \beta \right)}   % alpha/beta  for slp computation
+\newcommand{\triad}[6][]{\ensuremath{{}_{#2}^{#3}{\mathbb{#4}_{#1}}_{#5}^{\,#6}}}
+\newcommand{\triadd}[5]{\ensuremath{{}_{#1}^{#2}{\mathbb{#3}}_{#4}^{\,#5}}}
+\newcommand{\triadt}[5]{\ensuremath{{}_{#1}^{#2}{\tilde{\mathbb{#3}}}_{#4}^{\,#5}}}
+\newcommand{\rtriad}[2][]{\ensuremath{\triad[#1]{i}{k}{#2}{i_p}{k_p}}}
+\newcommand{\rtriadt}[1]{\ensuremath{\triadt{i}{k}{#1}{i_p}{k_p}}}
+\newcommand{\Alts}{{A}}
+\newcommand{\Alt}{{A^{lT}}}
+\newcommand{\rMLt}[1][i]{\tilde{r}_{\mathrm{ML}\,#1}}
+\newcommand{\rML}[1][i]{r_{\mathrm{ML}\,#1}}
+\newcommand{\mygstrut}[2]{\rule[#1 em]{0pt}{#2 em}}
+\newcommand{\mystrut}{\rule[-.9 em]{0pt}{1.79 em}}