Changeset 1224

Timestamp:

2008-11-26T14:52:28+01:00 (15 years ago)

Author:

gm

Message:

minor corrections in the Chapters from Steven + gm see ticket #283

Location:

trunk/DOC/TexFiles/Chapters

Files:

• Abstracts_Foreword.tex (modified) (2 diffs)
• Chap_DOM.tex (modified) (11 diffs)
• Chap_DYN.tex (modified) (11 diffs)
• Chap_LBC.tex (modified) (11 diffs)
• Chap_LDF.tex (modified) (9 diffs)
• Chap_MISC.tex (modified) (8 diffs)
• Chap_Model_Basics.tex (modified) (47 diffs)
• Chap_SBC.tex (modified) (4 diffs)
• Chap_TRA.tex (modified) (41 diffs)
• Chap_ZDF.tex (modified) (18 diffs)
• Introduction.tex (modified) (4 diffs)

trunk/DOC/TexFiles/Chapters/Abstracts_Foreword.tex

@@ -9 +9 @@
 
 \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 be a flexible tool for studying the ocean and its interactions with the others components of the earth climate system (atmosphere, sea-ice, biogeochemical tracers, ...) over a wide range of space and time scales. Prognostic variables are the three-dimensional velocity field, a linear or non-linear sea surface height, the temperature and the salinity. In the horizontal direction, the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step $z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a three-dimensional Arakawa C-type grid. Various physical choices are available to describe ocean physics, including TKE and KPP vertical physics. Within NEMO, the ocean is interfaced with a sea-ice model (LIM), passive tracer and biogeochemical models (TOP) and, via the OASIS coupler, with several atmospheric general circulation models. 
+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 
+be a flexible tool for studying the ocean and its interactions with the others components of 
+the earth climate system (atmosphere, sea-ice, biogeochemical tracers, ...) over a wide range 
+of space and time scales. Prognostic variables are the three-dimensional velocity field, a linear 
+or non-linear sea surface height, the temperature and the salinity. In the horizontal direction, 
+the model uses a curvilinear orthogonal grid and in the vertical direction, a full or partial step 
+$z$-coordinate, or $s$-coordinate, or a mixture of the two. The distribution of variables is a 
+three-dimensional Arakawa C-type grid. Various physical choices are available to describe 
+ocean physics, including TKE and KPP vertical physics. Within NEMO, the ocean is interfaced 
+with a sea-ice model (LIM v2 and v3), passive tracer and biogeochemical models (TOP) 
+and, via the OASIS coupler, with several atmospheric general circulation models. 
 
 % ================================================================
  \vspace{0.5cm}
 
-Le moteur oc\'{e}anique de NEMO (Nucleus for European Modelling of the Ocean) est un mod\`{e}le aux \'{e}quations primitives de la circulation oc\'{e}anique r\'{e}gionale et globale. Il se veut un outil flexible pour \'{e}tudier sur un vaste spectre spatiotemporel l'oc\'{e}an et ses interactions avec les autres composantes du syst\`{e}me climatique terrestre (atmosph\`{e}re, glace de mer, traceurs biog\'{e}ochimiques...). Les variables pronostiques sont le champ tridimensionnel de vitesse, une hauteur de la mer lin\'{e}aire ou non, la temperature et la salinit\'{e}. La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore une combinaison des deux. Diff\'{e}rents choix sont propos\'{e}s pour d\'{e}crire la physique oc\'{e}anique, incluant notamment des physiques verticales TKE et KPP . A travers l'infrastructure NEMO, l'oc\'{e}an est interfac\'{e} avec un mod\`{e}le de glace de mer, des mod\`{e}les biog\'{e}ochimiques et de traceur passif, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les de circulation g\'{e}n\'{e}rale atmosph\'{e}rique. 
+Le moteur oc\'{e}anique de NEMO (Nucleus for European Modelling of the Ocean) est un 
+mod\`{e}le aux \'{e}quations primitives de la circulation oc\'{e}anique r\'{e}gionale et globale. 
+Il se veut un outil flexible pour \'{e}tudier sur un vaste spectre spatiotemporel l'oc\'{e}an et ses 
+interactions avec les autres composantes du syst\`{e}me climatique terrestre (atmosph\`{e}re, 
+glace de mer, traceurs biog\'{e}ochimiques...). Les variables pronostiques sont le champ 
+tridimensionnel de vitesse, une hauteur de la mer lin\'{e}aire ou non, la temperature et la salinit\'{e}. 
+La distribution des variables se fait sur une grille C d'Arakawa tridimensionnelle utilisant une 
+coordonn\'{e}e verticale $z$ \`{a} niveaux entiers ou partiels, ou une coordonn\'{e}e s, ou encore 
+une combinaison des deux. Diff\'{e}rents choix sont propos\'{e}s pour d\'{e}crire la physique 
+oc\'{e}anique, incluant notamment des physiques verticales TKE et KPP. A travers l'infrastructure 
+NEMO, l'oc\'{e}an est interfac\'{e} avec un mod\`{e}le de glace de mer, des mod\`{e}les 
+biog\'{e}ochimiques et de traceur passif, et, via le coupleur OASIS, \`{a} plusieurs mod\`{e}les 
+de circulation g\'{e}n\'{e}rale atmosph\'{e}rique. 
 } 
 
@@ -32 +55 @@
 
  \vspace{0.5cm}
-Additional information can be found on www.locean-ipsl.upmc.fr/NEMO
+Additional information can be found on http://www.nemo-ocean.eu/
  \vspace{0.5cm}
 

trunk/DOC/TexFiles/Chapters/Chap_DOM.tex

@@ -73 +73 @@
 provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial 
 derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and 
-$\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. Discrete partial derivatives are formulated by the traditional, centred second order 
+$\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. 
+Discrete partial derivatives are formulated by the traditional, centred second order 
 finite difference approximation while the scale factors are chosen equal to their 
 local analytical value. An important point here is that the partial derivative of the 
@@ -262 +263 @@
 same $k$ index, in opposition to what is done in the horizontal plane where 
 it is the $T$-point and the nearest velocity points in the direction of the horizontal 
-axis that have the same $i$ or $j$ index (compare the dashed area in Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are chosen 
-to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} code 
-\emph{before all the vertical derivatives} of the discrete equations given in this 
-documentation.
+axis that have the same $i$ or $j$ index (compare the dashed area in 
+Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are 
+chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} 
+code \emph{before all the vertical derivatives} of the discrete equations given in 
+this documentation.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -601 +603 @@
 through statement functions, using parameters provided in the \textit{par\_oce.h90} file. 
 
-It is possible to define a simple regular vertical grid by giving zero stretching (\pp{ppacr=0}). In that case, the parameters \jp{jpk} (number of $w$-levels) and \pp{pphmax} (total ocean depth in meters) fully define the grid. 
+It is possible to define a simple regular vertical grid by giving zero stretching (\pp{ppacr=0}). 
+In that case, the parameters \jp{jpk} (number of $w$-levels) and \pp{pphmax} 
+(total ocean depth in meters) fully define the grid. 
 
 For climate-related studies it is often desirable to concentrate the vertical resolution 
-near the ocean surface. The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 
+near the ocean surface. The following function is proposed as a standard for a 
+$z$-coordinate (with either full or partial steps): 
 \begin{equation} \label{DOM_zgr_ana}
 \begin{split}
@@ -715 +720 @@
 one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been 
 defined in the absence of bathymetry. With partial steps, layers from 1 to 
-\jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1) is 
-allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the 
+\jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1) 
+is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the 
 maximum thickness allowed is $2*e_{3t}(jpk-1)$. This has to be kept in mind when 
 specifying the maximum depth \pp{pphmax} in partial steps: for example, with 
-\pp{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). Two 
-variables in the namdom namelist are used to define the partial step 
+\pp{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth 
+allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). 
+Two variables in the namdom namelist are used to define the partial step 
 vertical grid. The mimimum water thickness (in meters) allowed for a cell 
 partially filled with bathymetry at level jk is the minimum of \np{e3zpsmin} 
@@ -750 +756 @@
 surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean 
 depth, since a mixed step-like and bottom-following representation of the 
-topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided (\hf{zgr\_s} file) $h$ is a smooth envelope bathymetry and steps are used to represent sharp bathymetric gradients.
+topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided 
+(\hf{zgr\_s} file) $h$ is a smooth envelope bathymetry and steps are used to represent 
+sharp bathymetric gradients.
 
 A new flexible stretching function, modified from \citet{Song1994} is provided as an example:
@@ -763 +771 @@
 where $h_c$ is the thermocline depth and $\theta$ and $b$ are the surface and 
 bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 
-$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
+$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom 
+increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb] \label{Fig_sco_function}  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
-\caption{Examples of the stretching function applied to a sea mont; from left to right: surface, surface and bottom, and bottom intensified resolutions}
+\caption{Examples of the stretching function applied to a sea mont; from left to right: 
+surface, surface and bottom, and bottom intensified resolutions}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -859 +869 @@
 well as the implications in terms of starting or restarting a model 
 simulation. Note that the time stepping is generally performed in a one step 
-operation. With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in time for each term separately.
-%%%
-\gmcomment{ STEVEN  suggest separately instead of successively...  wrong?}
-%%%
+operation. With such a complex and nonlinear system of equations it would be 
+dangerous to let a prognostic variable evolve in time for each term separately.
 
 The three level scheme requires three arrays for each prognostic variables. 
@@ -896 +904 @@
 to diverge into a physical and a computational mode. Time splitting can 
 be controlled through the use of an Asselin time filter (first designed by 
-\citep{Robert1966} and more comprehensively studied by \citet{Asselin1972}), or by 
-periodically reinitialising the leapfrog solution through a single 
+\citep{Robert1966} and more comprehensively studied by \citet{Asselin1972}), 
+or by periodically reinitialising the leapfrog solution through a single 
 integration step with a two-level scheme. In \NEMO we follow the first 
 strategy:
@@ -996 +1004 @@
 \right.
 \end{equation}
-where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. The linear constraint \eqref{Eq_DOM_nxt_euler_stability} is a necessary condition, but not sufficient. If it is not satisfied, even mildly, then the model soon becomes wildly unstable. The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient.
+where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is 
+the mixing coefficient. The linear constraint \eqref{Eq_DOM_nxt_euler_stability} 
+is a necessary condition, but not sufficient. If it is not satisfied, even mildly, 
+then the model soon becomes wildly unstable. The instability can be removed 
+by either reducing the length of the time steps or reducing the mixing coefficient.
 
 For the vertical diffusion terms, a forward time differencing scheme can be 
@@ -1032 +1044 @@
 \right]
 \end{equation}
-where RHS is the right hand side of the equation except for the vertical diffusion term. We rewrite \eqref{Eq_DOM_nxt_imp} as:
+where RHS is the right hand side of the equation except for the vertical diffusion term. 
+We rewrite \eqref{Eq_DOM_nxt_imp} as:
 \begin{equation} \label{Eq_DOM_nxt_imp_mat}
 -c(k+1)\;u^{t+1}(k+1)+d(k)\;u^{t+1}(k)-\;c(k)\;u^{t+1}(k-1) \equiv b(k)
@@ -1075 +1088 @@
 gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be 
 added in the restart file to ensure an exact restartability. This is done only optionally 
-via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of restart file can be obtained when the restartability is not a key issue (operational oceanography or ensemble simulation for seasonal forcast).
+via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of 
+restart file can be obtained when the restartability is not a key issue (operational 
+oceanography or ensemble simulation for seasonal forcast).
 %%%
 \gmcomment{add here how to force the restart to contain only one time step for operational purposes}

trunk/DOC/TexFiles/Chapters/Chap_DYN.tex

@@ -8 +8 @@
 % add a figure for  dynvor ens, ene latices
 
-
+%\vspace{2.cm}
 $\ $\newline      %force an empty line
 
@@ -59 +59 @@
 MISC correspond to "extracting tendency terms" or "vorticity balance"?}
 
+$\ $\newline    % force a new ligne
 % ================================================================
 % Coriolis and Advection terms: vector invariant form
@@ -70 +71 @@
 The vector invariant form of the momentum equations is the one most 
 often used in applications of \NEMO ocean model. The flux form option 
-(see next section) has been recently introduced in version $2$. 
-Coriolis and momentum 
-advection terms are evaluated using a leapfrog scheme, $i.e.$ the velocity 
-appearing in these expressions is centred in time (\textit{now} velocity). 
+(see next section) has been introduced since version $2$. 
+Coriolis and momentum advection terms are evaluated using a leapfrog 
+scheme, $i.e.$ the velocity appearing in these expressions is centred in 
+time (\textit{now} velocity). 
 At the lateral boundaries either free slip, no slip or partial slip boundary 
 conditions are applied following Chap.\ref{LBC}.
@@ -88 +89 @@
 
 Different discretisations of the vorticity term (\textit{ln\_dynvor\_xxx}=.true.) are 
-available: conserving potential enstrophy of horizontally non-divergent flow; 
-conserving horizontal kinetic energy; or conserving potential enstrophy for the 
-relative vorticity term and horizontal kinetic energy for the planetary vorticity term 
-(see  Appendix~\ref{Apdx_C}). The vorticity terms are given below for the general 
-case, but note that in the full step $z$-coordinate (\key{zco} is defined), 
-$e_{3u} =e_{3v} =e_{3f}$ so that the vertical scale factors disappear.
+available: conserving potential enstrophy of horizontally non-divergent flow ; 
+conserving horizontal kinetic energy ; or conserving potential enstrophy for the 
+relative vorticity term and horizontal kinetic energy for the planetary vorticity 
+term (see  Appendix~\ref{Apdx_C}). The vorticity terms are given below for the 
+general case, but note that in the full step $z$-coordinate (\key{zco} is defined), 
+$e_{3u} =e_{3v} =e_{3f}$ so that the vertical scale factors disappear. They are 
+all computed in dedicated routines that can be found in the \mdl{dynvor} module.
 
 %-------------------------------------------------------------
@@ -108 +110 @@
 \left\{ 
 \begin{aligned}
-{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2}    \\
-{-\frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}  & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j}  
+{+\frac{1}{e_{1u} } } & {\overline {\left( { \frac{\zeta +f}{e_{3f} }} \right)} }^{\,i} 
+                                & {\overline{\overline {\left( {e_{1v} e_{3v} v} \right)}} }^{\,i, j+1/2}    \\
+{- \frac{1}{e_{2v} } } & {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)} }^{\,j}  
+                                & {\overline{\overline {\left( {e_{2u} e_{3u} u} \right)}} }^{\,i+1/2, j}  
 \end{aligned} 
  \right.
@@ -123 +127 @@
 kinetic energy but not the global enstrophy. It is given by:
 \begin{equation} \label{Eq_dynvor_ene}
-\left\{ {
-\begin{aligned}
-{+\frac{1}{e_{1u} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
-\;\overline {\left( {e_{1v} e_{3v} v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
-{-\frac{1}{e_{2v} }\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
-\;\overline {\left( {e_{2u} e_{3u} u} \right)} ^{\,j+1/2}} }^{\,i} }
-\end{aligned} 
-} \right.
+\left\{   \begin{aligned}
+{+\frac{1}{e_{1u}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
+                            \;  \overline {\left( {e_{1v} e_{3v} v} \right)} ^{\,i+1/2}} }^{\,j} }    \\
+{- \frac{1}{e_{2v}}\; {\overline {\left( {\frac{\zeta +f}{e_{3f} }} \right)
+                            \;  \overline {\left( {e_{2u} e_{3u} u} \right)} ^{\,j+1/2}} }^{\,i} }
+\end{aligned}    \right.
 \end{equation} 
 
@@ -230 +232 @@
 } \right.
 \end{equation} 
-where $a$, $b$, $c$ and $d$ are triad combinations of the neighbouring 
-potential vorticities (Fig. \ref{Fig_DYN_een_triad}): 
+where $a$, $b$, $c$ and $d$ are the following triad combinations of the 
+neighbouring potential vorticities (Fig.~\ref{Fig_DYN_een_triad}): 
 \begin{equation} \label{Eq_een_triads}
 \left\{ 
@@ -377 +379 @@
 
 The scheme is non diffusive (i.e. conserves the kinetic energy) but dispersive 
-($i.e.$ it may create false extrema). It is therefore notoriously noisy and must 
-be used in conjunction with an explicit diffusion operator to produce a sensible 
-solution. The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, so $u$ and $v$ are the \emph{now} 
-velocities.
+($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be 
+used in conjunction with an explicit diffusion operator to produce a sensible solution. 
+The associated time-stepping is performed using a leapfrog scheme in conjunction 
+with an Asselin time-filter, so $u$ and $v$ are the \emph{now} velocities.
 
 %-------------------------------------------------------------
@@ -716 +718 @@
 by $t-\Delta t$, $t, t+\Delta t$, and $t+2\Delta t$. 
 The curved line represents a leap-frog time step, and the smaller time 
-steps $N \Delta t_e=\frac{3}{2}\Delta t$ are denoted by the zig-zag line. The vertically 
-integrated forcing \textbf{M}(t) computed at the model time step $t$ 
+steps $N \Delta t_e=\frac{3}{2}\Delta t$ are denoted by the zig-zag line. 
+The vertically integrated forcing \textbf{M}(t) computed at the model time step $t$ 
 represents the interaction between the external and internal motions. 
-While keeping \textbf{M} and freshwater forcing field fixed, a 
-leap-frog integration carries the external mode variables (surface height and vertically integrated velocity) from $t$ to $t+\frac{3}{2} \Delta t$ using N external time steps of length $\Delta t_e$. 
-Time averaging the external fields over the $\frac{2}{3}N+1$ time steps (endpoints 
-included) centers the vertically integrated velocity and the sea surface height at the model timestep $t+\Delta t$. These averaged values are used to update \textbf{M}(t) with both the surface pressure gradient and the Coriolis force. 
-A baroclinic leap-frog time step carries the surface height to The model time stepping scheme can then be achieved by 
-$t+\Delta t$ using the convergence of the time averaged vertically integrated 
-velocity taken from baroclinic time step t. }
-%%%
-\gmcomment{STEVEN: what does convergence mean in this context?}
-%%%
+While keeping \textbf{M} and freshwater forcing field fixed, a leap-frog 
+integration carries the external mode variables (surface height and vertically 
+integrated velocity) from $t$ to $t+\frac{3}{2} \Delta t$ using N external time 
+steps of length $\Delta t_e$. Time averaging the external fields over the 
+$\frac{2}{3}N+1$ time steps (endpoints included) centers the vertically integrated 
+velocity and the sea surface height at the model timestep $t+\Delta t$. 
+These averaged values are used to update \textbf{M}(t) with both the surface 
+pressure gradient and the Coriolis force, therefore providing the $t+\Delta t$
+velocity.  The model time stepping scheme can then be achieved by a baroclinic 
+leap-frog time step that carries the surface height from $t-\Delta t$ to $t+\Delta t$.  }
 \end{center}
 \end{figure}
@@ -988 +990 @@
 
 The turbulent flux of momentum at the bottom of the ocean is specified through 
-a bottom friction parameterization (see \S\ref{ZDF_bfr})
+a bottom friction parameterisation (see \S\ref{ZDF_bfr})
 
 % ================================================================
@@ -1061 +1063 @@
 \end{equation} 
 
-Note that in the $z$-coordinate with full step (\key{zco} is defined), $e_{3u} =e_{3v} =e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}.
+Note that in the $z$-coordinate with full step (\key{zco} is defined), 
+$e_{3u} =e_{3v} =e_{3f}$ so that they cancel in \eqref{Eq_divcur_div}.
 
 Note also that whereas the vorticity have the same discrete expression in $z$- 

trunk/DOC/TexFiles/Chapters/Chap_LBC.tex

@@ -5 +5 @@
 \label{LBC}
 \minitoc
+
+\newpage
+$\ $\newline    % force a new ligne
+
 
 %gm% add here introduction to this chapter
@@ -146 +150 @@
 \label{LBC_jperio}
 
-At the model domain boundaries several choices are offered: closed, cyclic east-west, south symmetric across the equator, a north-fold, and combination closed-north fold or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh. 
+At the model domain boundaries several choices are offered: closed, cyclic east-west, 
+south symmetric across the equator, a north-fold, and combination closed-north fold 
+or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -198 +204 @@
 \label{LBC_north_fold}
 
-The north fold boundary condition has been introduced in order to handle the north boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere. \colorbox{yellow}{to be completed...}
+The north fold boundary condition has been introduced in order to handle the north 
+boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere. 
+\colorbox{yellow}{to be completed...}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_North_Fold_T}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_North_Fold_T.pdf}
-\caption {North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds to the inner domain mask (see text). }
+\caption {North fold boundary with a $T$-point pivot and cyclic east-west boundary condition 
+($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds to the inner 
+domain mask (see text). }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -240 +250 @@
 and the solving of the elliptic equation associated with the surface pressure 
 gradient computation (delocalization over the whole horizontal domain). 
-Therefore, a pencil strategy is used for the data sub-structuration \gmcomment{no 
-idea what this means!}: the 3D initial domain is laid out on local processor 
+Therefore, a pencil strategy is used for the data sub-structuration 
+\gmcomment{no idea what this means!}
+: the 3D initial domain is laid out on local processor 
 memories following a 2D horizontal topological splitting. Each sub-domain 
 computes its own surface and bottom boundary conditions and has a side 
@@ -249 +260 @@
 phase starts: each processor sends to its neighbouring processors the update 
 values of the points corresponding to the interior overlapping area to its 
-neighbouring sub-domain (i.e. the innermost of the two overlapping rows). The communication is done through message passing. Usually the parallel virtual 
+neighbouring sub-domain (i.e. the innermost of the two overlapping rows). 
+The communication is done through message passing. Usually the parallel virtual 
 language, PVM, is used as it is a standard language available on  nearly  all 
 MPP computers. More specific languages (i.e. computer dependant languages) 
@@ -272 +284 @@
  \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in 
  \mdl{par\_oce}). Each processor is independent and without message passing 
- or synchronous process \gmcomment{how does a synchronous process relate to this?}, 
+ or synchronous process 
+ \gmcomment{how does a synchronous process relate to this?}, 
  programs run alone and access just its own local memory. For this reason, the 
  main model dimensions are now the local dimensions of the subdomain (pencil) 
@@ -286 +299 @@
 where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis.
 
-\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and no east-west cyclic boundary conditions.}
+\colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and 
+no east-west cyclic boundary conditions.}
 
 One also defines variables nldi and nlei which correspond to the internal 
@@ -309 +323 @@
 \item       nbondi =  2    no splitting following the i-axis.
 \end{itemize}
-   During the simulation, processors exchange data with their neighbours. If there is effectively a neighbour, the processor receives variables from this processor on its overlapping row, and sends the data issued from internal domain corresponding to the overlapping row of the other processor.
+During the simulation, processors exchange data with their neighbours. 
+If there is effectively a neighbour, the processor receives variables from this 
+processor on its overlapping row, and sends the data issued from internal 
+domain corresponding to the overlapping row of the other processor.
        
-
 \colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos }
 
 
-
-The OPA model computes equation terms with the help of mask arrays (0 on land 
+The \NEMO model computes equation terms with the help of mask arrays (0 on land 
 points and 1 on sea points). It is easily readable and very efficient in the context of 
 a computer with vectorial architecture. However, in the case of a scalar processor, 
@@ -332 +347 @@
 nono, noea,...) so that the land-only processors are not taken into account. 
 
-\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp routine should be suppressed from the code.}
-
-When land processors are eliminated, the value corresponding to these locations in the model output files is zero. Note that this is a problem for a mesh output file written by such a model configuration, because model users often divide by the scale factors ($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be best not to eliminate land processors when running the model especially to write the mesh files as outputs (when \np{nmsh} namelist parameter differs from 0).
-\gmcomment{Steven : dont understand this, no land processor means no output file covering this part of globe; its only when files are stitched together into one that you can leave a hole}
+\colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp 
+routine should be suppressed from the code.}
+
+When land processors are eliminated, the value corresponding to these locations in 
+the model output files is zero. Note that this is a problem for a mesh output file written 
+by such a model configuration, because model users often divide by the scale factors 
+($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be 
+best not to eliminate land processors when running the model especially to write the 
+mesh files as outputs (when \np{nmsh} namelist parameter differs from 0).
+\gmcomment{Steven : dont understand this, no land processor means no output file 
+covering this part of globe; its only when files are stitched together into one that you 
+can leave a hole}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -383 +406 @@
 NEMO (originally in OPA8.2). It allows the user to 
 \begin{itemize}
-\item tell the model that a boundary is ''open'' and not closed by a wall, for example by modifying the calculation of the divergence of velocity there;
-\item impose values of tracers and velocities at that boundary (values which may be taken from a climatology): this is the``fixed OBC'' option. 
-\item calculate boundary values by a sophisticated algorithm combining radiation and relaxation (``radiative OBC'' option)
+\item tell the model that a boundary is ''open'' and not closed by a wall, for example 
+by modifying the calculation of the divergence of velocity there;
+\item impose values of tracers and velocities at that boundary (values which may 
+be taken from a climatology): this is the``fixed OBC'' option. 
+\item calculate boundary values by a sophisticated algorithm combining radiation 
+and relaxation (``radiative OBC'' option)
 \end{itemize}
 
@@ -529 +555 @@
 to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the 
 small domain can be created by using the following netcdf utility on the global files: 
-ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities, see http://nco.sourceforge.net). The open boundary files can be constructed using ncks 
+ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities, see http://nco.sourceforge.net). 
+The open boundary files can be constructed using ncks 
 commands, following table~\ref{Tab_obc_ind}. 
 

trunk/DOC/TexFiles/Chapters/Chap_LDF.tex

@@ -75 +75 @@
 
 \subsubsection{Horizontally Varying Mixing Coefficients (\key{ldftra\_c2d} and \key{ldfdyn\_c2d})}
-By default the horizontal variation of the eddy coefficient depends on the local mesh size and the type of operator used:
+By default the horizontal variation of the eddy coefficient depends on the local mesh 
+size and the type of operator used:
 \begin{equation} \label{Eq_title}
   A_l = \left\{     
@@ -84 +85 @@
 \quad \text{comments}
 \end{equation}
-where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked ocean domain, and $A_o^l$ is the \np{ahm0} (momentum) or \np{aht0} (tracer) namelist parameter. This variation is intended to reflect the lesser need for subgrid scale eddy mixing where the grid size is smaller in the domain. It was introduced in the context of the DYNAMO modelling project \citep{Willebrand2001}. 
+where $e_{max}$ is the maximum of $e_1$ and $e_2$ taken over the whole masked 
+ocean domain, and $A_o^l$ is the \np{ahm0} (momentum) or \np{aht0} (tracer) 
+namelist parameter. This variation is intended to reflect the lesser need for subgrid 
+scale eddy mixing where the grid size is smaller in the domain. It was introduced in 
+the context of the DYNAMO modelling project \citep{Willebrand2001}. 
 %%%
-\gmcomment { not only that! stability reasons: with non uniform grid size, it is common to face a blow up of the model due to to large diffusive coefficient compare to the smallest grid size... this is especially true for bilaplacian (to be added in the text!)  }
+\gmcomment { not only that! stability reasons: with non uniform grid size, it is common 
+to face a blow up of the model due to to large diffusive coefficient compare to the smallest 
+grid size... this is especially true for bilaplacian (to be added in the text!)  }
 
 Other formulations can be introduced by the user for a given configuration. 
@@ -150 +157 @@
 
 %%%
-\gmcomment{  we should emphasize here that the implementation is a rather old one. Better work can be achieved by using \citet{Griffies1998, Griffies2004} iso-neutral scheme. }
+\gmcomment{  we should emphasize here that the implementation is a rather old one. 
+Better work can be achieved by using \citet{Griffies1998, Griffies2004} iso-neutral scheme. }
 
 A direction for lateral mixing has to be defined when the desired operator does 
@@ -227 +235 @@
 \end{equation}
 
-%gm% rewrite this as the explanation in not very clear !!!
+%gm% rewrite this as the explanation is not very clear !!!
 %In practice, \eqref{Eq_ldfslp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \eqref{Eq_ldfslp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 
 
@@ -262 +270 @@
 
 %gm% 
-Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for the constraint on iso-neutral fluxes. Following \citet{Griffies2004}, instead of specifying directly that there is a zero neutral diffusive flux of locally referenced potential density, we stay in the $T$-$S$ plane and consider the balance between the neutral direction diffusive fluxes of potential temperature and salinity:
+Note: The solution for $s$-coordinate passes trough the use of different 
+(and better) expression for the constraint on iso-neutral fluxes. Following 
+\citet{Griffies2004}, instead of specifying directly that there is a zero neutral 
+diffusive flux of locally referenced potential density, we stay in the $T$-$S$ 
+plane and consider the balance between the neutral direction diffusive fluxes 
+of potential temperature and salinity:
 \begin{equation}
 \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S)
@@ -323 +336 @@
 the effect of an horizontal background mixing. 
 
-Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, contrary to the \citet{Griffies1998} operator which has that property. 
+Nevertheless, this iso-neutral operator does not ensure that variance cannot increase, 
+contrary to the \citet{Griffies1998} operator which has that property. 
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -332 +346 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-%There are three additional questions about the slope calculation. First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes. Second, numerical stability issues also require a bound on slopes. Third, the question of boundary condition specified on slopes...
+%There are three additional questions about the slope calculation. 
+%First the expression for the rotation tensor has been obtain assuming the "small slope" approximation, so a bound has to be imposed on slopes. 
+%Second, numerical stability issues also require a bound on slopes. 
+%Third, the question of boundary condition specified on slopes...
 
 %from griffies: chapter 13.1....
@@ -338 +355 @@
 
 
-In addition and also for numerical stability reasons \citep{Cox1987, Griffies2004}, 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, Griffies2004}, 
+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, Griffies2004}, the slopes must also 
@@ -404 +422 @@
 an eddy induced tracer advection term is added, 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} is used in $z$-coordinates, and the sum \eqref{Eq_ldfslp_geo}  
+mixing, the slopes used here are referenced to the geopotential surfaces, $i.e.$ 
+\eqref{Eq_ldfslp_geo} is used in $z$-coordinates, and the sum \eqref{Eq_ldfslp_geo}  
 + \eqref{Eq_ldfslp_iso} in $s$-coordinates. The eddy induced velocity is given by: 
 \begin{equation} \label{Eq_ldfeiv}

trunk/DOC/TexFiles/Chapters/Chap_MISC.tex

@@ -40 +40 @@
 \textit{Top}: using partially open cells. The meridional scale factor at $v$-point 
 is reduced on both sides of the strait to account for the real width of the strait 
-(about 20 km). Note that the scale factors of the strait $T$-point remains unchanged. \textit{Bottom}: using viscous boundary layers. The four fmask parameters 
+(about 20 km). Note that the scale factors of the strait $T$-point remains unchanged. 
+\textit{Bottom}: using viscous boundary layers. The four fmask parameters 
 along the strait coastlines are set to a value larger than 4, $i.e.$ "strong" no-slip 
 case (see Fig.\ref{Fig_LBC_shlat}) creating a large viscous boundary layer 
@@ -124 +125 @@
 \key{cfg\_1d} CPP key. This 1D model is a very useful tool  \textit{(a)} to learn 
 about the physics and numerical treatment of vertical mixing processes ; \textit{(b)} 
-to investigate suitable parameterizations of unresolved turbulence (wind steering, 
+to investigate suitable parameterisations of unresolved turbulence (wind steering, 
 langmuir circulation, skin layers) ; \textit{(c)} to compare the behaviour of different 
 vertical mixing schemes  ; \textit{(d)} to perform sensitivity studies on the vertical 
@@ -254 +255 @@
 size. This allows a very large model domain to be used, just by changing the domain 
 size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step or the physical 
-parameterizations. 
+parameterisations. 
 
 
@@ -429 +430 @@
 {\textbf{ x}}^{n+1}=\text{inf} _{{\textbf{ y}}={\textbf{ x}}^n+\alpha^n \,{\textbf{ d}}^n} \,\phi ({\textbf{ y}})\;\;\Leftrightarrow \;\;\frac{d\phi }{d\alpha}=0
 \end{equation*}
-and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the value that minimises the functional: 
+and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the 
+value that minimises the functional: 
 \begin{equation*}
 \alpha ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle  / \langle {\textbf{ A d}^n, \textbf{d}^n} \rangle
@@ -439 +441 @@
 product linked to \textbf{A}. Expressing the condition 
 $\langle \textbf{A d}^n, \textbf{d}^{n-1} \rangle = 0$ the value of $\beta ^n$ is found:
- $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle  / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. As a result, the errors $ \textbf{r}^n$ form an orthogonal 
+ $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle  / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$. 
+ As a result, the errors $ \textbf{r}^n$ form an orthogonal 
 base for the canonic dot product while the descent vectors $\textbf{d}^n$ form 
 an orthogonal base for the dot product linked to \textbf{A}. The resulting 
@@ -497 +500 @@
 The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the 
 canonical dot product the following one is used: 
-${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, 
-and if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ are substituted to \textbf{b} and \textbf{A} \citep{Madec1988}. 
-In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for \textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of 
+${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and 
+if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ 
+are substituted to \textbf{b} and \textbf{A} \citep{Madec1988}. 
+In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for 
+\textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of 
 \eqref{Eq_pmat} is in fact given by \eqref{Eq_solmat_p} and thus the matrix and 
 right hand side are computed independently from the solver used.
@@ -609 +614 @@
 %       Tracer/Dynamics Trends
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra}, \key{diatrddyn})}
+\subsection[Tracer/Dynamics Trends (\key{trdlmd}, \textbf{key\_diatrd...})]
+                  {Tracer/Dynamics Trends (\key{trdlmd}, \key{diatrdtra}, \key{diatrddyn})}
 \label{MISC_tratrd}
 
@@ -636 +642 @@
 
 The on-line computation of floats adevected either by the three dimensional velocity 
-field or constraint to remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. The algorithm used is based on 
+field or constraint to remain at a given depth ($w = 0$ in the computation) have been 
+introduced in the system during the CLIPPER project. The algorithm used is based on 
 the work of \cite{Blanke_Raynaud_JPO97}. (see also the web site describing the off-line 
 use of this marvellous diagnostic tool (http://stockage.univ-brest.fr/~grima/Ariane/).

trunk/DOC/TexFiles/Chapters/Chap_Model_Basics.tex

@@ -27 +27 @@
 velocity, plus the following additional assumptions made from scale considerations:
 
-\textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to be spheres so that gravity (local vertical) is parallel to the earth's radius
+\textit{(1) spherical earth approximation: }the geopotential surfaces are assumed to 
+be spheres so that gravity (local vertical) is parallel to the earth's radius
 
 \textit{(2) thin-shell approximation: }the ocean depth is neglected compared to the earth's radius
 
-\textit{(3) turbulent closure hypothesis: }the turbulent fluxes (which represent the effect of small scale processes on the large-scale) are expressed in terms of large-scale features
-
-\textit{(4) Boussinesq hypothesis:} density variations are neglected except in their contribution to the buoyancy force
-
-\textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a balance between the vertical pressure gradient and the buoyancy force (this removes convective processes from 
-the initial Navier-Stokes equations and so convective processes must be parameterized instead)
-
-\textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity vector is assumed to be zero.
-
-Because the gravitational force is so dominant in the equations of large-scale motions, it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ (the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane), $T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) vector system provides the following six equations (namely the momentum balance, the hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation equations and an equation of state):
+\textit{(3) turbulent closure hypothesis: }the turbulent fluxes (which represent the effect 
+of small scale processes on the large-scale) are expressed in terms of large-scale features
+
+\textit{(4) Boussinesq hypothesis:} density variations are neglected except in their 
+contribution to the buoyancy force
+
+\textit{(5) Hydrostatic hypothesis: }the vertical momentum equation is reduced to a 
+balance between the vertical pressure gradient and the buoyancy force (this removes 
+convective processes from the initial Navier-Stokes equations and so convective processes 
+must be parameterized instead)
+
+\textit{(6) Incompressibility hypothesis: }the three dimensional divergence of the velocity 
+vector is assumed to be zero.
+
+Because the gravitational force is so dominant in the equations of large-scale motions, 
+it is useful to choose an orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked 
+to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two 
+vectors orthogonal to \textbf{k}, $i.e.$ tangent to the geopotential surfaces. Let us define 
+the following variables: \textbf{U} the vector velocity, $\textbf{U}=\textbf{U}_h + w\, \textbf{k}$ 
+(the subscript $h$ denotes the local horizontal vector, $i.e.$ over the (\textbf{i},\textbf{j}) plane), 
+$T$ the potential temperature, $S$ the salinity, \textit{$\rho $} the \textit{in situ} density. 
+The vector invariant form of the primitive equations in the (\textbf{i},\textbf{j},\textbf{k}) 
+vector system provides the following six equations (namely the momentum balance, the 
+hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation 
+equations and an equation of state):
 \begin{subequations} \label{Eq_PE}
   \begin{equation}     \label{Eq_PE_dyn}
@@ -65 +81 @@
   \end{equation}
 \end{subequations}
-where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, $t$ is the time, $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by the equation of state (\ref{Eq_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, $f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 
-${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterizations of small-scale physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$ surface forcing terms. Their nature and formulation are discussed in \S\ref{PE_zdf_ldf} and page \S\ref{PE_boundary_condition}.
+where $\nabla$ is the generalised derivative vector operator in $(\bf i,\bf j, \bf k)$ directions, 
+$t$ is the time, $z$ is the vertical coordinate, $\rho $ is the \textit{in situ} density given by 
+the equation of state (\ref{Eq_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 
+$f=2 \bf \Omega \cdot \bf k$ is the Coriolis acceleration (where $\bf \Omega$ is the Earth's 
+angular velocity vector), and $g$ is the gravitational acceleration. 
+${\rm {\bf D}}^{\rm {\bf U}}$, $D^T$ and $D^S$ are the parameterisations of small-scale 
+physics for momentum, temperature and salinity, and ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ 
+and $F^S$ surface forcing terms. Their nature and formulation are discussed in 
+\S\ref{PE_zdf_ldf} and page \S\ref{PE_boundary_condition}.
 
 .
@@ -76 +99 @@
 \label{PE_boundary_condition}
 
-An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, the continental margins, the sea ice and the atmosphere. However, some of these fluxes are so weak that even on climatic time scales of thousands of years they can be neglected. In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and the other components of the earth system.
+An ocean is bounded by complex coastlines, bottom topography at its base and an air-sea 
+or ice-sea interface at its top. These boundaries can be defined by two surfaces, $z=-H(i,j)$ 
+and $z=\eta(i,j,k,t)$, where $H$ is the depth of the ocean bottom and $\eta$ is the height 
+of the sea surface. Both $H$ and $\eta$ are usually referenced to a given surface, $z=0$, 
+chosen as a mean sea surface (Fig.~\ref{Fig_ocean_bc}). Through these two boundaries, 
+the ocean can exchange fluxes of heat, fresh water, salt, and momentum with the solid earth, 
+the continental margins, the sea ice and the atmosphere. However, some of these fluxes are 
+so weak that even on climatic time scales of thousands of years they can be neglected. 
+In the following, we briefly review the fluxes exchanged at the interfaces between the ocean 
+and the other components of the earth system.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!ht] \label{Fig_ocean_bc}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_I_ocean_bc.pdf}
-\caption{The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. Both $H$ and $\eta $ are referenced to $z=0$.}
+\caption{The ocean is bounded by two surfaces, $z=-H(i,j)$ and $z=\eta(i,j,k,t)$, where $H$ 
+is the depth of the sea floor and $\eta$ the height of the sea surface. Both $H$ and $\eta $ 
+are referenced to $z=0$.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -87 +121 @@
 
 \begin{description}
-\item[Land - ocean interface:] the major flux between continental margins and the ocean is a mass exchange of fresh water through river runoff. Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. It can be neglected for short range integrations but has to be taken into account for long term integrations as it influences the characteristics of water masses formed (especially at high latitudes). It is required in order to close the water cycle of the climate system. It is usually specified as a fresh water flux at the air-sea interface in the vicinity of river mouths.
-\item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small, except in special areas of little extent. They are usually neglected in the model \footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling ($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world ocean (see \ref{TRA_bbc}).}. The boundary condition is thus set to no flux of heat and salt across solid boundaries. For momentum, the situation is different. There is no flow across solid boundaries, $i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, the bottom velocity is parallel to solid boundaries). This kinematic boundary condition can be expressed as:
+\item[Land - ocean interface:] the major flux between continental margins and the ocean is 
+a mass exchange of fresh water through river runoff. Such an exchange modifies the sea 
+surface salinity especially in the vicinity of major river mouths. It can be neglected for short 
+range integrations but has to be taken into account for long term integrations as it influences 
+the characteristics of water masses formed (especially at high latitudes). It is required in order 
+to close the water cycle of the climate system. It is usually specified as a fresh water flux at 
+the air-sea interface in the vicinity of river mouths.
+\item[Solid earth - ocean interface:] heat and salt fluxes through the sea floor are small, 
+except in special areas of little extent. They are usually neglected in the model 
+\footnote{In fact, it has been shown that the heat flux associated with the solid Earth cooling 
+($i.e.$the geothermal heating) is not negligible for the thermohaline circulation of the world 
+ocean (see \ref{TRA_bbc}).}. 
+The boundary condition is thus set to no flux of heat and salt across solid boundaries. 
+For momentum, the situation is different. There is no flow across solid boundaries, 
+$i.e.$ the velocity normal to the ocean bottom and coastlines is zero (in other words, 
+the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 
+can be expressed as:
 \begin{equation} \label{Eq_PE_w_bbc}
 w = -{\rm {\bf U}}_h \cdot  \nabla _h \left( H \right)
 \end{equation}
-In addition, the ocean exchanges momentum with the earth through frictional processes. Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification depends on the nature of the physical parameterization used for ${\rm {\bf D}}^{\rm {\bf U}}$ in \eqref{Eq_PE_dyn}. It is discussed in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9.
-\item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux of fresh water PE  (the precipitation minus evaporation budget) leads to: 
+In addition, the ocean exchanges momentum with the earth through frictional processes. 
+Such momentum transfer occurs at small scales in a boundary layer. It must be parameterized 
+in terms of turbulent fluxes using bottom and/or lateral boundary conditions. Its specification 
+depends on the nature of the physical parameterisation used for ${\rm {\bf D}}^{\rm {\bf U}}$ 
+in \eqref{Eq_PE_dyn}. It is discussed in \S\ref{PE_zdf}, page~\pageref{PE_zdf}.% and Chap. III.6 to 9.
+\item[Atmosphere - ocean interface:] the kinematic surface condition plus the mass flux 
+of fresh water PE  (the precipitation minus evaporation budget) leads to: 
 \begin{equation} \label{Eq_PE_w_sbc}
 w = \frac{\partial \eta }{\partial t} 
@@ -99 +153 @@
     + P-E
 \end{equation}
-The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) leads to the continuity of pressure across the interface $z=\eta$. The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat.
-\item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water and momentum. The sea surface temperature is constrained to be at the freezing point at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected.
+The dynamic boundary condition, neglecting the surface tension (which removes capillary 
+waves from the system) leads to the continuity of pressure across the interface $z=\eta$. 
+The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat.
+\item[Sea ice - ocean interface:] the ocean and sea ice exchange heat, salt, fresh water 
+and momentum. The sea surface temperature is constrained to be at the freezing point 
+at the interface. Sea ice salinity is very low ($\sim4-6 \,psu$) compared to those of the 
+ocean ($\sim34 \,psu$). The cycle of freezing/melting is associated with fresh water and 
+salt fluxes that cannot be neglected.
 \end{description}
 
@@ -116 +176 @@
 \label{PE_p_formulation}
 
-The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: $p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}), assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}). The hydrostatic pressure is then given by:
+The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at a 
+reference geopotential surface ($z=0$) and a hydrostatic pressure $p_h$ such that: 
+$p(i,j,k,t)=p_s(i,j,t)+p_h(i,j,k,t)$. The latter is computed by integrating (\ref{Eq_PE_hydrostatic}), 
+assuming that pressure in decibars can be approximated by depth in meters in (\ref{Eq_PE_eos}). 
+The hydrostatic pressure is then given by:
 \begin{equation} \label{Eq_PE_pressure}
 p_h \left( {i,j,z,t} \right)
- = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,z} \right)\;d\varsigma } 
-\end{equation}
- Two strategies can be considered for the surface pressure term: $(a)$ introduce of a new variable $\eta$, the free-surface elevation, for which a prognostic equation can be established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, one solution of the free-surface elevation consists of the excitation of external gravity waves. The flow is barotropic and the surface moves up and down with gravity as the restoring force. The phase speed of such waves is high (some hundreds of metres per second) so that the time step would have to be very short if they were present in the model. The latter strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$ the sea surface is the surface $z=0$. This well known approximation increases the surface wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic Rossby or planetary waves). In the present release of \NEMO, both strategies are still available. They are further described in the next two sub-sections.
+ = \int_{\varsigma =z}^{\varsigma =0} {g\;\rho \left( {T,S,\varsigma} \right)\;d\varsigma } 
+\end{equation}
+ Two strategies can be considered for the surface pressure term: $(a)$ introduce of a 
+ new variable $\eta$, the free-surface elevation, for which a prognostic equation can be 
+ established and solved; $(b)$ assume that the ocean surface is a rigid lid, on which the 
+ pressure (or its horizontal gradient) can be diagnosed. When the former strategy is used, 
+ one solution of the free-surface elevation consists of the excitation of external gravity waves. 
+ The flow is barotropic and the surface moves up and down with gravity as the restoring force. 
+ The phase speed of such waves is high (some hundreds of metres per second) so that 
+ the time step would have to be very short if they were present in the model. The latter 
+ strategy filters out these waves since the rigid lid approximation implies $\eta=0$, $i.e.$ 
+ the sea surface is the surface $z=0$. This well known approximation increases the surface 
+ wave speed to infinity and modifies certain other longwave dynamics ($e.g.$ barotropic 
+ Rossby or planetary waves). In the present release of \NEMO, both strategies are still available. 
+ They are further described in the next two sub-sections.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -129 +205 @@
 \label{PE_free_surface}
 
-In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced which describes the shape of the air-sea interface. This variable is solution of a prognostic equation which is established by forming the vertical average of the kinematic surface condition (\ref{Eq_PE_w_bbc}):
+In the free surface formulation, a variable $\eta$, the sea-surface height, is introduced 
+which describes the shape of the air-sea interface. This variable is solution of a 
+prognostic equation which is established by forming the vertical average of the kinematic 
+surface condition (\ref{Eq_PE_w_bbc}):
 \begin{equation} \label{Eq_PE_ssh}
 \frac{\partial \eta }{\partial t}=-D+P-E
@@ -137 +216 @@
 and using (\ref{Eq_PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$.
 
-Allowing the air-sea interface to move introduces the external gravity waves (EGWs) as a class of solution of the primitive equations. These waves are barotropic because of hydrostatic assumption, and their phase speed is quite high. Their time scale is short with respect to the other processes described by the primitive equations.
-
-Three choices can be made regarding the implementation of the free surface in the model, depending on the physical processes of interest. 
+Allowing the air-sea interface to move introduces the external gravity waves (EGWs) 
+as a class of solution of the primitive equations. These waves are barotropic because 
+of hydrostatic assumption, and their phase speed is quite high. Their time scale is 
+short with respect to the other processes described by the primitive equations.
+
+Three choices can be made regarding the implementation of the free surface in the model, 
+depending on the physical processes of interest. 
 
 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction 
-with the baroclinic structure of the ocean (internal waves) possibly in 
-shallow seas, then a non linear free surface is the most appropriate. This 
-means that no approximation is made in (\ref{Eq_PE_ssh}) and that the variation of the 
-ocean volume is fully taken into account. Note that in order to study the 
-fast time scales associated with EGWs it is necessary to minimize time 
-filtering effects (use an explicit time scheme with very small time step, or 
-a split-explicit scheme with reasonably small time step, see \S\ref{DYN_spg_exp} or
-\S\ref{DYN_spg_ts}.
+with the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 
+then a non linear free surface is the most appropriate. This means that no 
+approximation is made in (\ref{Eq_PE_ssh}) and that the variation of the ocean 
+volume is fully taken into account. Note that in order to study the fast time scales 
+associated with EGWs it is necessary to minimize time filtering effects (use an 
+explicit time scheme with very small time step, or a split-explicit scheme with 
+reasonably small time step, see \S\ref{DYN_spg_exp} or \S\ref{DYN_spg_ts}.
 
 $\bullet$ If one is not interested in EGW but rather sees them as high frequency 
@@ -163 +245 @@
 external waves are removed from the system. 
 
-The filtering of EGWs in models with a free surface is usually a matter of discretisation of the temporal derivatives, using the time splitting method \citep{Killworth1991, Zhang1992} or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach developed by \citet{Roullet2000}: the damping of EGWs is ensured by introducing an additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes: 
+The filtering of EGWs in models with a free surface is usually a matter of discretisation 
+of the temporal derivatives, using the time splitting method \citep{Killworth1991, Zhang1992} 
+or the implicit scheme \citep{Dukowicz1994}. In \NEMO, we use a slightly different approach 
+developed by \citet{Roullet2000}: the damping of EGWs is ensured by introducing an 
+additional force in the momentum equation. \eqref{Eq_PE_dyn} becomes: 
 \begin{equation} \label{Eq_PE_flt}
 \frac{\partial {\rm {\bf U}}_h }{\partial t}= {\rm {\bf M}}
@@ -169 +255 @@
 - g \ T_c \nabla \left( \widetilde{\rho} \ \partial_t \eta \right) 
 \end{equation}
-where $T_c$, is a parameter with dimensions of time which characterizes the force, $\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, non-linear and viscous terms in \eqref{Eq_PE_dyn}.
-
-The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate, $i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than $T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs can be damped by choosing $T_c > \Delta t$. \citet{Roullet2000} demonstrate that (\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which has to be computed implicitly. This is not surprising since the use of a large time step has a necessarily numerical cost. Two gains arise in comparison with the previous formulations. Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as soon as $T_c > \Delta t$.
-
-When the variations of free surface elevation are small compared to the thickness of the first model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized by \citet{Roullet2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the conservation of salt in the model. With the nonlinear free surface equation, the time evolution of the total salt content is 
+where $T_c$, is a parameter with dimensions of time which characterizes the force, 
+$\widetilde{\rho} = \rho / \rho_o$ is the dimensionless density, and $\rm {\bf M}$ 
+represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
+non-linear and viscous terms in \eqref{Eq_PE_dyn}.
+
+The new force can be interpreted as a diffusion of vertically integrated volume flux divergence. 
+The time evolution of $D$ is thus governed by a balance of two terms, $-g$ \textbf{A} $\eta$ 
+and $g \, T_c \,$ \textbf{A} $D$, associated with a propagative regime and a diffusive regime 
+in the temporal spectrum, respectively. In the diffusive regime, the EGWs no longer propagate, 
+$i.e.$ they are stationary and damped. The diffusion regime applies to the modes shorter than 
+$T_c$. For longer ones, the diffusion term vanishes. Hence, the temporally unresolved EGWs 
+can be damped by choosing $T_c > \Delta t$. \citet{Roullet2000} demonstrate that 
+(\ref{Eq_PE_flt}) can be integrated with a leap frog scheme except the additional term which 
+has to be computed implicitly. This is not surprising since the use of a large time step has a 
+necessarily numerical cost. Two gains arise in comparison with the previous formulations. 
+Firstly, the damping of EGWs can be quantified through the magnitude of the additional term. 
+Secondly, the numerical scheme does not need any tuning. Numerical stability is ensured as 
+soon as $T_c > \Delta t$.
+
+When the variations of free surface elevation are small compared to the thickness of the first 
+model layer, the free surface equation (\ref{Eq_PE_ssh}) can be linearized. As emphasized 
+by \citet{Roullet2000} the linearization of (\ref{Eq_PE_ssh}) has consequences on the 
+conservation of salt in the model. With the nonlinear free surface equation, the time evolution 
+of the total salt content is 
 \begin{equation} \label{Eq_PE_salt_content}
-\frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} =\int\limits_S 
-{S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds}
-\end{equation}
-where $S$ is the salinity, and the total salt is integrated over the whole ocean volume $D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) is satisfied, so that the salt is perfectly conserved. When the free surface equation is linearized, \citet{Roullet2000} show that the total salt content integrated in the fixed volume $D$ (bounded by the surface $z=0$) is no longer conserved:
+    \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} 
+                        =\int\limits_S {S\;(-\frac{\partial \eta }{\partial t}-D+P-E)\;ds}
+\end{equation}
+where $S$ is the salinity, and the total salt is integrated over the whole ocean volume 
+$D_\eta$ bounded by the time-dependent free surface. The right hand side (which is an 
+integral over the free surface) vanishes when the nonlinear equation (\ref{Eq_PE_ssh}) 
+is satisfied, so that the salt is perfectly conserved. When the free surface equation is 
+linearized, \citet{Roullet2000} show that the total salt content integrated in the fixed 
+volume $D$ (bounded by the surface $z=0$) is no longer conserved:
 \begin{equation} \label{Eq_PE_salt_content_linear}
-\frac{\partial }{\partial t}\int\limits_D {S\;dv} =-\int\limits_S 
-{S\;\frac{\partial \eta }{\partial t}ds} 
-\end{equation}
-
-The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions \citep{Roullet2000}. It can be significant when the freshwater forcing is not balanced and the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} results in a decrease of the salinity in the fixed volume $D$. Even in that case though, the total salt integrated in the variable volume $D_{\eta}$ varies much less, since (\ref{Eq_PE_salt_content_linear}) can be rewritten as 
+         \frac{\partial }{\partial t}\int\limits_D {S\;dv} 
+               = - \int\limits_S {S\;\frac{\partial \eta }{\partial t}ds} 
+\end{equation}
+
+The right hand side of (\ref{Eq_PE_salt_content_linear}) is small in equilibrium solutions 
+\citep{Roullet2000}. It can be significant when the freshwater forcing is not balanced and 
+the globally averaged free surface is drifting. An increase in sea surface height \textit{$\eta $} 
+results in a decrease of the salinity in the fixed volume $D$. Even in that case though, 
+the total salt integrated in the variable volume $D_{\eta}$ varies much less, since 
+(\ref{Eq_PE_salt_content_linear}) can be rewritten as 
 \begin{equation} \label{Eq_PE_salt_content_corrected}
 \frac{\partial }{\partial t}\int\limits_{D\eta } {S\;dv} 
@@ -191 +306 @@
 \end{equation}
 
-Although the total salt content is not exactly conserved with the linearized free surface, its variations are driven by correlations of the time variation of surface salinity with the sea surface height, which is a negligible term. This situation contrasts with 
-the case of the rigid lid approximation (following section) in which case freshwater forcing is represented by a virtual salt flux, leading to a spurious source of salt at the ocean surface \citep{Roullet2000}.
+Although the total salt content is not exactly conserved with the linearized free surface, 
+its variations are driven by correlations of the time variation of surface salinity with the 
+sea surface height, which is a negligible term. This situation contrasts with the case of 
+the rigid lid approximation (following section) in which case freshwater forcing is 
+represented by a virtual salt flux, leading to a spurious source of salt at the ocean 
+surface \citep{Roullet2000}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -200 +319 @@
 \label{PE_rigid_lid}
 
-With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface is equal to zero. From the continuity equation \eqref{Eq_PE_continuity} and the kinematic condition at the bottom \eqref{Eq_PE_w_bbc} (no flux across the bottom), it can be shown that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the overbar indicates a vertical average over the whole water column, i.e. from $z=-H$, the ocean bottom, to $z=0$ , the rigid-lid). Thus, $\rm {\bf \bar {U}}_h$ can be derived from a volume transport streamfunction $\psi$:
+With the rigid lid approximation, we assume that the ocean surface ($z=0$) is a rigid lid 
+on which a pressure $p_s$ is exerted. This implies that the vertical velocity at the surface 
+is equal to zero. From the continuity equation \eqref{Eq_PE_continuity} and the kinematic 
+condition at the bottom \eqref{Eq_PE_w_bbc} (no flux across the bottom), it can be shown 
+that the vertically integrated flow $H{\rm {\bf \bar {U}}}_h$ is non-divergent (where the 
+overbar indicates a vertical average over the whole water column, i.e. from $z=-H$, 
+the ocean bottom, to $z=0$ , the rigid-lid). Thus, $\rm {\bf \bar {U}}_h$ can be derived 
+from a volume transport streamfunction $\psi$:
 \begin{equation} \label{Eq_PE_u_psi}
 \overline{\vect{U}}_h =\frac{1}{H}\left(   \vect{k} \times \nabla \psi   \right)
 \end{equation}
 
-As $p_s$ does not depend on depth, its horizontal gradient is obtained by forming the vertical average of \eqref{Eq_PE_dyn} and using \eqref{Eq_PE_u_psi}:
+As $p_s$ does not depend on depth, its horizontal gradient is obtained by forming the 
+vertical average of \eqref{Eq_PE_dyn} and using \eqref{Eq_PE_u_psi}:
 
 \begin{equation} \label{Eq_PE_u_barotrope}
@@ -214 +341 @@
 \end{equation}
 
-Here ${\rm {\bf M}} = \left( M_u,M_v \right)$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, nonlinear and viscous terms in \eqref{Eq_PE_dyn}. The time derivative of $\psi $ is the solution of an elliptic equation which is obtained from the vertical component of the curl of (\ref{Eq_PE_u_barotrope}):
+Here ${\rm {\bf M}} = \left( M_u,M_v \right)$ represents the collected contributions of the 
+Coriolis, hydrostatic pressure gradient, nonlinear and viscous terms in \eqref{Eq_PE_dyn}. 
+The time derivative of $\psi $ is the solution of an elliptic equation which is obtained from 
+the vertical component of the curl of (\ref{Eq_PE_u_barotrope}):
 \begin{equation} \label{Eq_PE_psi}
 \left[   {\nabla \times \left[ {\frac{1}{H} \vect{\bf k} 
@@ -221 +351 @@
 \end{equation}
 
-Using the proper boundary conditions, \eqref{Eq_PE_psi} can be solved to find $\partial_t \psi$ and thus using \eqref{Eq_PE_u_barotrope} the horizontal surface pressure gradient. It should be noted that $p_s$ can be computed by taking the divergence of \eqref{Eq_PE_u_barotrope} and solving the resulting elliptic equation. Thus the surface pressure is a diagnostic quantity that can be recovered for analysis purposes.
-
-A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. The boundary condition on velocity is that there is no flow normal to a solid wall, $i.e.$ the coastlines are streamlines. Therefore \eqref{Eq_PE_psi} is solved with the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each coastline of the same continent or of the same island. When all the coastlines are connected (there are no islands), the constant value of $\partial_t \psi$ along the coast can be arbitrarily chosen to be zero. When islands are present in the domain, the value of the barotropic streamfunction will generally be different for each island and for the continent, and will vary with respect to time. So the boundary condition is: $\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), where $Q$ is the number of islands present in the domain and $\mu_n$ is a time dependent variable. A time evolution equation of the unknown $\mu_n$ can be found by evaluating the circulation of the time derivative of the vertical average (barotropic) velocity field along a closed contour around each island. Since the circulation of a gradient field along a closed contour is zero, from \eqref{Eq_PE_u_barotrope} we have:
+Using the proper boundary conditions, \eqref{Eq_PE_psi} can be solved to find $\partial_t \psi$ 
+and thus using \eqref{Eq_PE_u_barotrope} the horizontal surface pressure gradient. 
+It should be noted that $p_s$ can be computed by taking the divergence of 
+\eqref{Eq_PE_u_barotrope} and solving the resulting elliptic equation. Thus the surface 
+pressure is a diagnostic quantity that can be recovered for analysis purposes.
+
+A difficulty lies in the determination of the boundary condition on $\partial_t \psi$. 
+The boundary condition on velocity is that there is no flow normal to a solid wall, 
+$i.e.$ the coastlines are streamlines. Therefore \eqref{Eq_PE_psi} is solved with 
+the following Dirichlet boundary condition: $\partial_t \psi$ is constant along each 
+coastline of the same continent or of the same island. When all the coastlines are 
+connected (there are no islands), the constant value of $\partial_t \psi$ along the 
+coast can be arbitrarily chosen to be zero. When islands are present in the domain, 
+the value of the barotropic streamfunction will generally be different for each island 
+and for the continent, and will vary with respect to time. So the boundary condition is: 
+$\psi=0$ along the continent and $\psi=\mu_n$ along island $n$ ($1 \leq n \leq Q$), 
+where $Q$ is the number of islands present in the domain and $\mu_n$ is a time 
+dependent variable. A time evolution equation of the unknown $\mu_n$ can be found 
+by evaluating the circulation of the time derivative of the vertical average (barotropic) 
+velocity field along a closed contour around each island. Since the circulation of a 
+gradient field along a closed contour is zero, from \eqref{Eq_PE_u_barotrope} we have:
 \begin{equation} \label{Eq_PE_isl_circulation}
 \oint_n {\frac{1}{H}\left[ {{\rm {\bf k}}\times \nabla \left( 
@@ -256 +404 @@
 \right)_{1\leqslant n\leqslant Q} ={\rm {\bf B}}
 \end{equation}
-where \textbf{A} is a $Q  \times Q$ matrix and \textbf{B} is a time dependent vector. As \textbf{A} is independent of time, it can be calculated and inverted once. The time derivative of the streamfunction when islands are present is thus given by:
+where \textbf{A} is a $Q  \times Q$ matrix and \textbf{B} is a time dependent vector. 
+As \textbf{A} is independent of time, it can be calculated and inverted once. The time 
+derivative of the streamfunction when islands are present is thus given by:
 \begin{equation} \label{Eq_PE_psi_isl_dt}
 \frac{\partial \psi }{\partial t}=\frac{\partial \psi _o }{\partial 
@@ -277 +427 @@
 \label{PE_tensorial}
 
-In many ocean circulation problems, the flow field has regions of enhanced dynamics ($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). The representation of such dynamical processes can be improved by specifically increasing the model resolution in these regions. As well, it may be convenient to use a lateral boundary-following coordinate system to better represent coastal dynamics. Moreover, the common geographical coordinate system has a singular point at the North Pole that cannot be easily treated in a global model without filtering. A solution consists of introducing an appropriate coordinate transformation that shifts the singular point onto land \citep{MadecImb1996, Murray1996}. As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems. An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. This formalism is suited to any multidimensional curvilinear coordinate system. Ocean modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth approximation), with preservation of the local vertical. Here we give the simplified equations for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey of the conservation laws of fluid dynamics.
-
-Let (\textbf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (Fig.\ref{Fig_referential}). Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea level (Fig.\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, the three scale factors:
+In many ocean circulation problems, the flow field has regions of enhanced dynamics 
+($i.e.$ surface layers, western boundary currents, equatorial currents, or ocean fronts). 
+The representation of such dynamical processes can be improved by specifically increasing 
+the model resolution in these regions. As well, it may be convenient to use a lateral 
+boundary-following coordinate system to better represent coastal dynamics. Moreover, 
+the common geographical coordinate system has a singular point at the North Pole that 
+cannot be easily treated in a global model without filtering. A solution consists of introducing 
+an appropriate coordinate transformation that shifts the singular point onto land 
+\citep{MadecImb1996, Murray1996}. As a consequence, it is important to solve the primitive 
+equations in various curvilinear coordinate systems. An efficient way of introducing an 
+appropriate coordinate transform can be found when using a tensorial formalism. 
+This formalism is suited to any multidimensional curvilinear coordinate system. Ocean 
+modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth 
+approximation), with preservation of the local vertical. Here we give the simplified equations 
+for this particular case. The general case is detailed by \citet{Eiseman1980} in their survey 
+of the conservation laws of fluid dynamics.
+
+Let (\textbf{i},\textbf{j},\textbf{k}) be a set of orthogonal curvilinear coordinates on the sphere 
+associated with the positively oriented orthogonal set of unit vectors (\textbf{i},\textbf{j},\textbf{k}) 
+linked to the earth such that \textbf{k} is the local upward vector and (\textbf{i},\textbf{j}) are 
+two vectors orthogonal to \textbf{k}, $i.e.$ along geopotential surfaces (Fig.\ref{Fig_referential}). 
+Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined 
+by the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and the distance from the centre of 
+the earth $a+z(k)$ where $a$ is the earth's radius and $z$ the altitude above a reference sea 
+level (Fig.\ref{Fig_referential}). The local deformation of the curvilinear coordinate system is 
+given by $e_1$, $e_2$ and $e_3$, the three scale factors:
 \begin{equation} \label{Eq_scale_factors}
 \begin{aligned}
@@ -295 +468 @@
 \begin{figure}[!tb] \label{Fig_referential}  \begin{center}
 \includegraphics[width=0.60\textwidth]{./TexFiles/Figures/Fig_I_earth_referential.pdf}
-\caption{the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear coordinate system (\textbf{i},\textbf{j},\textbf{k}). }
+\caption{the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
+coordinate system (\textbf{i},\textbf{j},\textbf{k}). }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by $a$ in (\ref{Eq_scale_factors}) (thin-shell approximation). The resulting horizontal scale factors $e_1$, $e_2$  are independent of $k$ while the vertical scale factor is a single function of $k$ as \textbf{k} is parallel to \textbf{z}. The scalar and vector operators that appear in the primitive equations (Eqs. \eqref{Eq_PE_dyn} to \eqref{Eq_PE_eos}) can be written in the tensorial form, invariant in any orthogonal horizontal curvilinear coordinate system transformation:
+Since the ocean depth is far smaller than the earth's radius, $a+z$, can be replaced by 
+$a$ in (\ref{Eq_scale_factors}) (thin-shell approximation). The resulting horizontal scale 
+factors $e_1$, $e_2$  are independent of $k$ while the vertical scale factor is a single 
+function of $k$ as \textbf{k} is parallel to \textbf{z}. The scalar and vector operators that 
+appear in the primitive equations (Eqs. \eqref{Eq_PE_dyn} to \eqref{Eq_PE_eos}) can 
+be written in the tensorial form, invariant in any orthogonal horizontal curvilinear coordinate 
+system transformation:
 \begin{subequations} \label{Eq_PE_discrete_operators}
 \begin{equation} \label{Eq_PE_grad}
@@ -341 +521 @@
 \label{PE_zco_Eq}
 
-In order to express the Primitive Equations in tensorial formalism, it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using \eqref{Eq_PE_grad}) to \eqref{Eq_PE_lap_vector}. Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by:
+In order to express the Primitive Equations in tensorial formalism, it is necessary to compute 
+the horizontal component of the non-linear and viscous terms of the equation using 
+\eqref{Eq_PE_grad}) to \eqref{Eq_PE_lap_vector}. 
+Let us set $\vect U=(u,v,w)={\vect{U}}_h +w\;\vect{k}$, the velocity in the $(i,j,k)$ coordinate 
+system and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity 
+field $\chi$, by:
 \begin{equation} \label{Eq_PE_curl_Uh}
 \zeta =\frac{1}{e_1 e_2 }\left[ {\frac{\partial \left( {e_2 \,v} 
@@ -353 +538 @@
 \end{equation}
 
-Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that $e_3$  is a function of the single variable $k$, the nonlinear term of \eqref{Eq_PE_dyn} can be transformed as follows:
+Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ 
+and that $e_3$  is a function of the single variable $k$, the nonlinear term of 
+\eqref{Eq_PE_dyn} can be transformed as follows:
 \begin{flalign*}
 &\left[ {\left( { \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
@@ -391 +578 @@
 \end{flalign*}
 
-The last term of the right hand side is obviously zero, and thus the nonlinear term of \eqref{Eq_PE_dyn} is written in the $(i,j,k)$ coordinate system:
+The last term of the right hand side is obviously zero, and thus the nonlinear term of 
+\eqref{Eq_PE_dyn} is written in the $(i,j,k)$ coordinate system:
 \begin{equation} \label{Eq_PE_vector_form}
 \left[ {\left( {  \nabla \times {\rm {\bf U}}    } \right) \times {\rm {\bf U}}
@@ -401 +589 @@
 \end{equation}
 
-This is the so-called \textit{vector invariant form} of the momentum advection term. For some purposes, it can be advantageous to write this term in the so-called flux form, $i.e.$ to write it as the divergence of fluxes. For example, the first component of \eqref{Eq_PE_vector_form} (the $i$-component) is transformed as follows:
+This is the so-called \textit{vector invariant form} of the momentum advection term. 
+For some purposes, it can be advantageous to write this term in the so-called flux form, 
+$i.e.$ to write it as the divergence of fluxes. For example, the first component of 
+\eqref{Eq_PE_vector_form} (the $i$-component) is transformed as follows:
 \begin{flalign*}
 &{ \begin{array}{*{20}l}
@@ -486 +677 @@
 \end{multline}
 
-The flux form has two terms, the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) and the second one is due to the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} term and can be viewed as a modification of the Coriolis parameter: 
+The flux form has two terms, the first one is expressed as the divergence of momentum 
+fluxes (hence the flux form name given to this formulation) and the second one is due to 
+the curvilinear nature of the coordinate system used. The latter is called the \emph{metric} 
+term and can be viewed as a modification of the Coriolis parameter: 
 \begin{equation} \label{Eq_PE_cor+metric}
 f \to f + \frac{1}{e_1 \; e_2}   \left(    v \frac{\partial e_2}{\partial i}
@@ -492 +686 @@
 \end{equation}
 
-Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$ and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of the Coriolis parameter $f \to f+(u/a) \tan \varphi$.
+Note that in the case of geographical coordinate, $i.e.$ when $(i,j) \to (\lambda ,\varphi )$ 
+and $(e_1 ,e_2) \to (a \,\cos \varphi ,a)$, we recover the commonly used modification of 
+the Coriolis parameter $f \to f+(u/a) \tan \varphi$.
 
 To sum up, the equations solved by the ocean model can be written in the following tensorial formalism:
@@ -545 +741 @@
 \end{multline}
 \end{subequations}
-where $\zeta$ is given by \eqref{Eq_PE_curl_Uh} and the surface pressure gradient formulation depends on the one of the free surface:
+where $\zeta$ is given by \eqref{Eq_PE_curl_Uh} and the surface pressure gradient formulation 
+depends on the one of the free surface:
 
 $*$ free surface formulation
@@ -566 +763 @@
 \end{equation}
 where ${\vect{M}}= \left( M_u,M_v \right)$ represents the collected contributions of nonlinear, 
-viscosity and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation:
+viscosity and hydrostatic pressure gradient terms in \eqref{Eq_PE_dyn_vect} and the overbar 
+indicates a vertical average over the whole water column ($i.e.$ from $z=-H$, the ocean bottom, 
+to $z=0$, the rigid-lid), and where the time derivative of $\psi$ is the solution of an elliptic equation:
 \begin{multline} \label{Eq_psi_total}
   \frac{\partial }{\partial i}\left[ {\frac{e_2 }{H\,e_1}\frac{\partial}{\partial i}
@@ -605 +804 @@
 \end{equation}
 
-The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid 
-scale parameterization used. It will be defined in \S\ref{PE_zdf}. The nature and formulation of ${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed in Chapter~\ref{SBC}.
+The expression of \textbf{D}$^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale 
+parameterisation used. It will be defined in \S\ref{PE_zdf}. The nature and formulation of 
+${\rm {\bf F}}^{\rm {\bf U}}$, $F^T$ and $F^S$, the surface forcing terms, are discussed 
+in Chapter~\ref{SBC}.
 
 \newpage 
@@ -617 +818 @@
 \begin{figure}[!b] \label{Fig_z_zstar}  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zstar.pdf}
-\caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate \citep{Adcroft_Campin_OM04} ).}
+\caption{(a) $z$-coordinate in linear free-surface case ; (b) $z-$coordinate in non-linear 
+free surface case (c) re-scaled height coordinate (become popular as the \textit{z*-}coordinate 
+\citep{Adcroft_Campin_OM04} ).}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 
-In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. These coordinates systems is presented in a report 
-\citep{Levier2007} available on the \NEMO web site. 
+In that case, the free surface equation is nonlinear, and the variations of volume are fully 
+taken into account. These coordinates systems is presented in a report \citep{Levier2007} 
+available on the \NEMO web site. 
 
 \gmcomment{
-The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation which allows one to deal with large amplitude free-surface
+The \textit{z*} coordinate approach is an unapproximated, non-linear free surface implementation 
+which allows one to deal with large amplitude free-surface
 variations relative to the vertical resolution \citep{Adcroft_Campin_OM04}. In
 the  \textit{z*} formulation, the variation of the column thickness due to sea-surface
@@ -642 +847 @@
 
 Since the vertical displacement of the free surface is incorporated in the vertical
-coordinate  \textit{z*}, the upper and lower boundaries are at fixed  \textit{z*} position,  $\textit{z*} = 0$ and  $\textit{z*} = ?H$ respectively. Also the divergence of the flow field is no longer zero as shown by the continuity equation:
+coordinate  \textit{z*}, the upper and lower boundaries are at fixed  \textit{z*} position,  
+$\textit{z*} = 0$ and  $\textit{z*} = ?H$ respectively. Also the divergence of the flow field 
+is no longer zero as shown by the continuity equation:
 
 $\frac{\partial r}{\partial t} = \nabla_{\textit{z*}} \cdot \left( r \; \rm{\bf U}_h \right)
@@ -662 +869 @@
 \subsection{Introduction}
 
-Several important aspects of the ocean circulation are influenced by bottom topography. Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, sills and channels that strongly constrain the path of water masses, but more subtle effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. Topographic Rossby waves can be excited and can interact with the mean current. In the $z-$coordinate system presented in the previous section (\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom and to large localized depth gradients associated with large localized vertical velocities. The response to such a velocity field often leads to numerical dispersion effects. One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) 
-
-The $s$-coordinate avoids the discretisation error in the depth field since the layers of computation are gradually adjusted with depth to the ocean bottom. Relatively small topographic features as well as  gentle, large-scale slopes of the sea floor in the deep ocean, which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, can easily be represented (with relatively low vertical resolution). A terrain-following model (hereafter $s-$model) also facilitates the modelling of the boundary layer flows over a large depth range, which in the framework of the $z$-model would require high vertical resolution over the whole depth range. Moreover, with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as the only boundaries of the domain (nomore lateral boundary condition to specify). Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a homogeneous ocean, it has strong limitations as soon as stratification is introduced. The main two problems come from the truncation error in the horizontal pressure gradient and a possibly increased diapycnal diffusion. The horizontal pressure force in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}),
+Several important aspects of the ocean circulation are influenced by bottom topography. 
+Of course, the most important is that bottom topography determines deep ocean sub-basins, 
+barriers, sills and channels that strongly constrain the path of water masses, but more subtle 
+effects exist. For example, the topographic $\beta$-effect is usually larger than the planetary 
+one along continental slopes. Topographic Rossby waves can be excited and can interact 
+with the mean current. In the $z-$coordinate system presented in the previous section 
+(\S\ref{PE_zco}), $z-$surfaces are geopotential surfaces. The bottom topography is 
+discretised by steps. This often leads to a misrepresentation of a gradually sloping bottom 
+and to large localized depth gradients associated with large localized vertical velocities. 
+The response to such a velocity field often leads to numerical dispersion effects. 
+One solution to strongly reduce this error is to use a partial step representation of bottom 
+topography instead of a full step one \cite{Pacanowski_Gnanadesikan_MWR98}. 
+Another solution is to introduce a terrain-following coordinate system (hereafter $s-$coordinate) 
+
+The $s$-coordinate avoids the discretisation error in the depth field since the layers of 
+computation are gradually adjusted with depth to the ocean bottom. Relatively small 
+topographic features as well as  gentle, large-scale slopes of the sea floor in the deep 
+ocean, which would be ignored in typical $z$-model applications with the largest grid 
+spacing at greatest depths, can easily be represented (with relatively low vertical resolution). 
+A terrain-following model (hereafter $s-$model) also facilitates the modelling of the 
+boundary layer flows over a large depth range, which in the framework of the $z$-model 
+would require high vertical resolution over the whole depth range. Moreover, with a 
+$s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface 
+as the only boundaries of the domain (nomore lateral boundary condition to specify). 
+Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a 
+homogeneous ocean, it has strong limitations as soon as stratification is introduced. 
+The main two problems come from the truncation error in the horizontal pressure 
+gradient and a possibly increased diapycnal diffusion. The horizontal pressure force 
+in $s$-coordinate consists of two terms (see Appendix~\ref{Apdx_A}),
 
 \begin{equation} \label{Eq_PE_p_sco}
@@ -671 +904 @@
 \end{equation}
 
-The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface and introduces a truncation error that is not present in a $z$-model. In the special case of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), \citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude of this truncation error. It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, and the finite difference scheme. This error limits the possible topographic slopes that a model can handle at a given horizontal and vertical resolution. This is a severe restriction for large-scale applications using realistic bottom topography. The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive. This problem can be at least partially overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec1996}. However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for a realistic bottom topography: a envelope topography is defined in $s$-coordinate on which a full or partial step bottom topography is then applied in order to adjust the model depth to the observed one (see \S\ref{DOM_zgr}.
-
-For numerical reasons a minimum of diffusion is required along the coordinate surfaces of any finite difference model. It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as well as for a $s$-model. 
-However, density varies more strongly on $s-$surfaces than on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. For example, imagine an isolated bump of topography in an ocean at rest with a horizontally uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column (i.e. the main thermocline) \citep{Madec1996}. An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}). 
-
-The $s-$coordinates introduced here \citep{Lott1990,Madec1996} differ mainly in two aspects from similar models:  it allows  a representation of bottom topography with mixed full or partial step-like/terrain following topography ; It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate.
+The second term in \eqref{Eq_PE_p_sco} depends on the tilt of the coordinate surface 
+and introduces a truncation error that is not present in a $z$-model. In the special case 
+of a $\sigma-$coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 
+\citet{Haney1991} and \citet{Beckmann1993} have given estimates of the magnitude 
+of this truncation error. It depends on topographic slope, stratification, horizontal and 
+vertical resolution, the equation of state, and the finite difference scheme. This error 
+limits the possible topographic slopes that a model can handle at a given horizontal 
+and vertical resolution. This is a severe restriction for large-scale applications using 
+realistic bottom topography. The large-scale slopes require high horizontal resolution, 
+and the computational cost becomes prohibitive. This problem can be at least partially 
+overcome by mixing $s$-coordinate and step-like representation of bottom topography \citep{Gerdes1993a,Gerdes1993b,Madec1996}. However, the definition of the model 
+domain vertical coordinate becomes then a non-trivial thing for a realistic bottom 
+topography: a envelope topography is defined in $s$-coordinate on which a full or 
+partial step bottom topography is then applied in order to adjust the model depth to 
+the observed one (see \S\ref{DOM_zgr}.
+
+For numerical reasons a minimum of diffusion is required along the coordinate surfaces 
+of any finite difference model. It causes spurious diapycnal mixing when coordinate 
+surfaces do not coincide with isoneutral surfaces. This is the case for a $z$-model as 
+well as for a $s$-model. However, density varies more strongly on $s-$surfaces than 
+on horizontal surfaces in regions of large topographic slopes, implying larger diapycnal 
+diffusion in a $s$-model than in a $z$-model. Whereas such a diapycnal diffusion in a 
+$z$-model tends to weaken horizontal density (pressure) gradients and thus the horizontal 
+circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 
+For example, imagine an isolated bump of topography in an ocean at rest with a horizontally 
+uniform stratification. Spurious diffusion along $s$-surfaces will induce a bump of isoneutral 
+surfaces over the topography, and thus will generate there a baroclinic eddy. In contrast, 
+the ocean will stay at rest in a $z$-model. As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below the strongly stratified portion of the water column 
+($i.e.$ the main thermocline) \citep{Madec1996}. An alternate solution consists of rotating 
+the lateral diffusive tensor to geopotential or to isoneutral surfaces (see \S\ref{PE_ldf}. 
+Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 
+strongly exceeding the stability limit of such a operator when it is discretized (see Chapter~\ref{LDF}). 
+
+The $s-$coordinates introduced here \citep{Lott1990,Madec1996} differ mainly in two 
+aspects from similar models:  it allows  a representation of bottom topography with mixed 
+full or partial step-like/terrain following topography ; It also offers a completely general 
+transformation, $s=s(i,j,z)$ for the vertical coordinate.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -683 +947 @@
 \subsection{The \textit{s-}coordinate Formulation}
 
-Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$ and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes $z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and $s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed equations is given in Appendix~\ref{Apdx_A}. Let us define the vertical scale factor by $e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the (\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by :
+Starting from the set of equations established in \S\ref{PE_zco} for the special case $k=z$ 
+and thus $e_3=1$, we introduce an arbitrary vertical coordinate $s=s(i,j,k,t)$, which includes 
+$z$-, \textit{z*}- and $\sigma-$coordinates as special cases ($s=z$, $s=\textit{z*}$, and 
+$s=\sigma=z/H$ or $=z/\left(H+\eta \right)$, resp.). A formal derivation of the transformed 
+equations is given in Appendix~\ref{Apdx_A}. Let us define the vertical scale factor by 
+$e_3=\partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), and the slopes in the 
+(\textbf{i},\textbf{j}) directions between $s-$ and $z-$surfaces by :
 \begin{equation} \label{Eq_PE_sco_slope}
 \sigma _1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s 
@@ -689 +959 @@
 \sigma _2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s
 \end{equation}
-We also introduce  $\omega $, a dia-surface velocity component, defined as the velocity relative to the moving $s$-surfaces and normal to them:
+We also introduce  $\omega $, a dia-surface velocity component, defined as the velocity 
+relative to the moving $s$-surfaces and normal to them:
 \begin{equation} \label{Eq_PE_sco_w}
 \omega  = w - e_3 \, \frac{\partial s}{\partial t} - \sigma _1 \,u - \sigma _2 \,v    \\
@@ -716 +987 @@
    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} \quad
 \end{multline}
-where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic pressure have the same expressions as in $z$-coordinates although they do not represent exactly the same quantities. $\omega$ is provided by the continuity equation (see Appendix~\ref{Apdx_A}):
+where the relative vorticity, \textit{$\zeta $}, the surface pressure gradient, and the hydrostatic 
+pressure have the same expressions as in $z$-coordinates although they do not represent 
+exactly the same quantities. $\omega$ is provided by the continuity equation 
+(see Appendix~\ref{Apdx_A}):
 
 \begin{equation} \label{Eq_PE_sco_continuity}
@@ -742 +1016 @@
 \end{multline}
 
-The equation of state has the same expression as in $z$-coordinate, and similar expressions are used for mixing and forcing terms.
+The equation of state has the same expression as in $z$-coordinate, and similar expressions 
+are used for mixing and forcing terms.
 
 \gmcomment{
@@ -768 +1043 @@
 It is usually called the subgrid scale physics. It must be emphasized that 
 this is the weakest part of the primitive equations, but also one of the 
-most important for long-term simulations as small scale processes \textit{in fine} balance 
-the surface input of kinetic energy and heat.
+most important for long-term simulations as small scale processes \textit{in fine} 
+balance the surface input of kinetic energy and heat.
 
 The control exerted by gravity on the flow induces a strong anisotropy 
-between the lateral and vertical motions. Therefore subgrid-scale physics  \textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in \eqref{Eq_PE_dyn}, \eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part  \textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part  \textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms and their underlying physics are briefly discussed in the next two subsections.
+between the lateral and vertical motions. Therefore subgrid-scale physics  
+\textbf{D}$^{\vect{U}}$, $D^{S}$ and $D^{T}$  in \eqref{Eq_PE_dyn}, 
+\eqref{Eq_PE_tra_T} and \eqref{Eq_PE_tra_S} are divided into a lateral part  
+\textbf{D}$^{l \vect{U}}$, $D^{lS}$ and $D^{lT}$ and a vertical part  
+\textbf{D}$^{vU}$, $D^{vS}$ and $D^{vT}$. The formulation of these terms 
+and their underlying physics are briefly discussed in the next two subsections.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -785 +1065 @@
 partially, but always parameterized. The vertical turbulent fluxes are 
 assumed to depend linearly on the gradients of large-scale quantities (for 
-example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, where $A^{vT}$ is an eddy coefficient). This formulation is 
+example, the turbulent heat flux is given by $\overline{T'w'}=-A^{vT} \partial_z \overline T$, 
+where $A^{vT}$ is an eddy coefficient). This formulation is 
 analogous to that of molecular diffusion and dissipation. This is quite 
 clearly a necessary compromise: considering only the molecular viscosity 
@@ -794 +1075 @@
 \begin{equation} \label{Eq_PE_zdf}
    \begin{split}
-{\vect{D}}^{v \vect{U}}
-&=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\         D^{vT} &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ ,
+{\vect{D}}^{v \vect{U}} &=\frac{\partial }{\partial z}\left( {A^{vm}\frac{\partial {\vect{U}}_h }{\partial z}} \right) \ , \\         
+D^{vT}                         &= \frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial T}{\partial z}} \right) \ ,
 \quad
 D^{vS}=\frac{\partial }{\partial z}\left( {A^{vT}\frac{\partial S}{\partial z}} \right)
    \end{split}
 \end{equation}
-where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified (see Chap.~\ref{SBC} and \ref{ZDF} and \S\ref{TRA_bbl}). All the vertical physics is embedded in the specification of the eddy coefficients. They can be assumed to be either constant, or function of the local fluid properties ($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a turbulent closure model. The choices available in \NEMO are discussed in \S\ref{ZDF}).
+where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, 
+respectively. At the sea surface and at the bottom, turbulent fluxes of momentum, heat 
+and salt must be specified (see Chap.~\ref{SBC} and \ref{ZDF} and \S\ref{TRA_bbl}). 
+All the vertical physics is embedded in the specification of the eddy coefficients. 
+They can be assumed to be either constant, or function of the local fluid properties 
+($e.g.$ Richardson number, Brunt-Vais\"{a}l\"{a} frequency...), or computed from a 
+turbulent closure model. The choices available in \NEMO are discussed in \S\ref{ZDF}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -821 +1108 @@
 lateral diffusive and dissipative operators are of second order. 
 Observations show that lateral mixing induced by mesoscale turbulence tends 
-to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987}) rather than across them. 
+to be along isopycnal surfaces (or more precisely neutral surfaces \cite{McDougall1987}) 
+rather than across them. 
 As the slope of neutral surfaces is small in the ocean, a common 
 approximation is to assume that the `lateral' direction is the horizontal, 
@@ -834 +1122 @@
 energy whereas potential energy is a main source of turbulence (through 
 baroclinic instabilities). \citet{Gent1990} have proposed a 
-parameterization of mesoscale eddy-induced turbulence which associates an 
+parameterisation of mesoscale eddy-induced turbulence which associates an 
 eddy-induced velocity to the isoneutral diffusion. Its mean effect is to 
 reduce the mean potential energy of the ocean. This leads to a formulation 
@@ -850 +1138 @@
 the model while not interfering with the solved 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
+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
 advection scheme. 
 
-All these parameterizations of subgrid scale physics present advantages and 
+All these parameterisations of subgrid scale physics present advantages and 
 drawbacks. There are not all available in \NEMO. In the $z$-coordinate 
 formulation, five options are offered for active tracers (temperature and 
 salinity): second order geopotential operator, second order isoneutral 
-operator, \citet{Gent1990} parameterization, fourth order 
-geopotential operator, and various slightly diffusive advection schemes. The same options are available for momentum, except 
-\citet{Gent1990} parameterization which only involves tracers. In the
+operator, \citet{Gent1990} parameterisation, fourth order 
+geopotential operator, and various slightly diffusive advection schemes. 
+The same options are available for momentum, except 
+\citet{Gent1990} parameterisation which only involves tracers. In the
 $s$-coordinate formulation, additional options are offered for tracers: second 
 order operator acting along $s-$surfaces, and for momentum: fourth order 
@@ -881 +1171 @@
 rotation between geopotential and $s$-surfaces, while it is only an approximation 
 for the rotation between isoneutral and $z$- or $s$-surfaces. Indeed, in the latter 
-case, two assumptions are made to simplify  \eqref{Eq_PE_iso_tensor} \citep{Cox1987}. First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity. Second, the two isoneutral directions of diffusion are assumed to be independent since the slopes are generally less than $10^{-2}$ in the ocean (see Appendix~\ref{Apdx_B}).
+case, two assumptions are made to simplify  \eqref{Eq_PE_iso_tensor} \citep{Cox1987}. 
+First, the horizontal contribution of the dianeutral mixing is neglected since the ratio 
+between iso and dia-neutral diffusive coefficients is known to be several orders of 
+magnitude smaller than unity. Second, the two isoneutral directions of diffusion are 
+assumed to be independent since the slopes are generally less than $10^{-2}$ in the 
+ocean (see Appendix~\ref{Apdx_B}).
 
 For \textit{geopotential} diffusion, $r_1$ and $r_2 $ are the slopes between the 
-geopotential and computational surfaces: in $z$-coordinates they are zero ($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ).
-
-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:
+geopotential and computational surfaces: in $z$-coordinates they are zero 
+($r_1 = r_2 = 0$) while in $s$-coordinate (including $\textit{z*}$ case) they are 
+equal to $\sigma _1$ and $\sigma _2$, respectively (see \eqref{Eq_PE_sco_slope} ).
+
+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:
 \begin{equation} \label{Eq_PE_iso_slopes}
 r_1 =\frac{e_3 }{e_1 }  \left( {\frac{\partial \rho }{\partial i}} \right)
@@ -894 +1192 @@
 \end{equation}
 
-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:
+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}
 D^{lT}=\nabla \cdot \left( {A^{lT}\;\Re \;\nabla T} \right)
            +\nabla \cdot \left( {{\vect{U}}^\ast \,T} \right)
 \end{equation}
-where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, eddy-induced transport velocity. This velocity field is defined by:
+where ${\vect{U}}^\ast =\left( {u^\ast ,v^\ast ,w^\ast } \right)$ is a non-divergent, 
+eddy-induced transport velocity. This velocity field is defined by:
 \begin{equation} \label{Eq_PE_eiv}
    \begin{split}
@@ -909 +1209 @@
    \end{split}
 \end{equation}
-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: 
+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: 
 \begin{align} \label{Eq_PE_slopes_eiv}
 \tilde{r}_n = \begin{cases}
@@ -918 +1221 @@
 \end{align}
 
-The normal component of the eddy induced velocity is zero at all the boundaries. this can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of the boundaries. The latter strategy is used in \NEMO (cf. Chap.~\ref{LDF}).
+The normal component of the eddy induced velocity is zero at all the boundaries. 
+This can be achieved in a model by tapering either the eddy coefficient or the slopes 
+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}
@@ -928 +1233 @@
  \end{equation}
 
-It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with the eddy diffusion coefficient correctly placed. 
+It is the second order operator given by \eqref{Eq_PE_iso_tensor} applied twice with 
+the eddy diffusion coefficient correctly placed. 
 
 
@@ -952 +1258 @@
 horizontal divergence fields (see Appendix~\ref{Apdx_C}). Unfortunately, it is not 
 available for geopotential diffusion in $s-$coordinates and for isoneutral 
-diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required). In these two cases, the $u$ and $v-$fields are considered as independent scalar fields, so that the diffusive operator is given by:
+diffusion in both $z$- and $s$-coordinates ($i.e.$ when a rotation is required). 
+In these two cases, the $u$ and $v-$fields are considered as independent scalar 
+fields, so that the diffusive operator is given by:
 \begin{equation} \label{Eq_PE_lapU_iso}
 \begin{split}
@@ -959 +1267 @@
  \end{split}
  \end{equation}
-where $\Re$ is given by  \eqref{Eq_PE_iso_tensor}. It is the same expression as those used for diffusive operator on tracers. It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or $\beta-$plane, not on the sphere. It is also a very good approximation in vicinity of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}.
+where $\Re$ is given by  \eqref{Eq_PE_iso_tensor}. It is the same expression as 
+those used for diffusive operator on tracers. It must be emphasised that such a 
+formulation is only exact in a Cartesian coordinate system, $i.e.$ on a $f-$ or 
+$\beta-$plane, not on the sphere. It is also a very good approximation in vicinity 
+of the Equator in a geographical coordinate system \citep{Lengaigne_al_JGR03}.
 
 \subsubsection{lateral fourth order momentum diffusive operator}
 
-As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} with the eddy viscosity coefficient correctly placed:
+As for tracers, the fourth order momentum diffusive operator along $z$ or $s$-surfaces 
+is a re-entering second order operator \eqref{Eq_PE_lapU} or \eqref{Eq_PE_lapU} 
+with the eddy viscosity coefficient correctly placed:
 
 geopotential diffusion in $z$-coordinate:

trunk/DOC/TexFiles/Chapters/Chap_SBC.tex

@@ -404 +404 @@
 %        Handling of ice-covered area
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Handling of ice-covered area}
+\subsection{Handling of ice-covered area  (\textit{sbcice\_...})}
 \label{SBC_ice-cover}
 
@@ -411 +411 @@
 depending on the value of the \np{nn{\_}ice} namelist parameter.  
 \begin{description}
-\item[nn{\_}ice = 0]  there will never be sea-ice in the computational domain. This is a typical namelist value used for tropical ocean domain. The surface fluxes are simply specified for an ice-free ocean. No specific things are done for sea-ice.
-\item[nn{\_}ice = 1]  sea-ice can exist in the computational domain, but no sea-ice model is used. An observed ice covered area is read in a file. Below this area, the SST is restored to the freezing point and the heat fluxes are set to $-4~W/m^2$ ($-2~W/m^2$) in the northern (southern) hemisphere. The associated modification of the freshwater fluxes are done in such a way that the change in buoyancy fluxes remains zero. This prevents deep convection to occur when trying to reach the freezing point (and so ice covered area condition) while the SSS is too large. This manner of managing sea-ice area, just by using si IF case, is usually referred as the \textit{ice-if} model. It can be found in the \mdl{sbcice{\_}if} module.
-\item[nn{\_}ice = 2 or more]  A full sea ice model is used. This model computes the ice-ocean fluxes, that are combined with the air-sea fluxes using the ice fraction of each model cell to provide the surface ocean fluxes. Note that the activation of a sea-ice model is is done by defining a CPP key (\key{lim2} or \key{lim3}). The activation automatically ovewrite the read value of nn{\_}ice to its appropriate value ($i.e.$ $2$ for LIM-2 and $3$ for LIM-3).
+\item[nn{\_}ice = 0]  there will never be sea-ice in the computational domain. 
+This is a typical namelist value used for tropical ocean domain. The surface fluxes 
+are simply specified for an ice-free ocean. No specific things is done for sea-ice.
+\item[nn{\_}ice = 1]  sea-ice can exist in the computational domain, but no sea-ice model 
+is used. An observed ice covered area is read in a file. Below this area, the SST is 
+restored to the freezing point and the heat fluxes are set to $-4~W/m^2$ ($-2~W/m^2$) 
+in the northern (southern) hemisphere. The associated modification of the freshwater 
+fluxes are done in such a way that the change in buoyancy fluxes remains zero. 
+This prevents deep convection to occur when trying to reach the freezing point 
+(and so ice covered area condition) while the SSS is too large. This manner of 
+managing sea-ice area, just by using si IF case, is usually referred as the \textit{ice-if} 
+model. It can be found in the \mdl{sbcice{\_}if} module.
+\item[nn{\_}ice = 2 or more]  A full sea ice model is used. This model computes the 
+ice-ocean fluxes, that are combined with the air-sea fluxes using the ice fraction of 
+each model cell to provide the surface ocean fluxes. Note that the activation of a 
+sea-ice model is is done by defining a CPP key (\key{lim2} or \key{lim3}). 
+The activation automatically ovewrite the read value of nn{\_}ice to its appropriate 
+value ($i.e.$ $2$ for LIM-2 and $3$ for LIM-3).
 \end{description}
 
@@ -428 +443 @@
 %-------------------------------------------------------------------------------------------------------------
 
+The river runoffs 
+
 It is convenient to introduce the river runoff in the model as a surface 
 fresh water flux. 
 
+
+%Griffies:  River runoff generally enters the ocean at a nonzero depth rather than through the surface. Many global models, however, have traditionally inserted river runoff to the top model cell. Such can become problematic numerically and physically when the top grid cells are reÞned to levels common in coastal modelling. Hence, more applications are now considering the input of runoff throughout a nonzero depth. Likewise, sea ice can melt at depth, thus necessitating a mass transport to occur within the ocean between the liquid and solid water masses.
+
 \colorbox{yellow}{Nevertheless, Pb of vertical resolution and increase of Kz in vicinity of }
 
 \colorbox{yellow}{river mouths{\ldots}}
 
-Control of the mean sea level
+%IF( ln_rnf ) THEN                                     ! increase diffusivity at rivers mouths
+%        DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + rn_avt_rnf * rnfmsk(:,:)   ;   END DO
+%ENDIF
+
+
 
 % -------------------------------------------------------------------------------------------------------------
@@ -444 +468 @@
 \label{SBC_fwb}
 
-To be written later...
-
-\gmcomment{The descrition of the technique used to control the freshwater budget has to be added here}
-
-
-
-
+For global ocean simulation it can be useful to introduce a control of the mean sea 
+level in order to prevent unrealistic drift of the sea surface height due to inaccuracy 
+in the freshwater fluxes. In \NEMO, two way of controlling the the freshwater budget. 
+\begin{description}
+\item[\np{nn\_fwb}=0]  no control at all. The mean sea level is free to drift, and will 
+certainly do so.
+\item[\np{nn\_fwb}=1]  global mean EMP set to zero at each model time step. 
+%Note that with a sea-ice model, this technique only control the mean sea level with linear free surface (\key{vvl} not defined) and no mass flux between ocean and ice (as it is implemented in the current ice-ocean coupling). 
+\item[\np{nn\_fwb}=2]  freshwater budget is adjusted from the previous year annual 
+mean budget which is read in the \textit{EMPave\_old.dat} file. As the model uses the 
+Boussinesq approximation, the annual mean fresh water budget is simply evaluated 
+from the change in the mean sea level at January the first and saved in the 
+\textit{EMPav.dat} file. 
+\end{description}
+
+% Griffies doc:
+% When running ocean-ice simulations, we are not explicitly representing land processes, such as rivers, catchment areas, snow accumulation, etc. However, to reduce model drift, it is important to balance the hydrological cycle in ocean-ice models. We thus need to prescribe some form of global normalization to the precipitation minus evaporation plus river runoff. The result of the normalization should be a global integrated zero net water input to the ocean-ice system over a chosen time scale. 
+%How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step, so that there is always a zero net input of water to the ocean-ice system. Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used to alter the subsequent yearÕs water budget in an attempt to damp the annual water imbalance. Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing. 
+%When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean and ice models when aiming to balance the hydrological cycle. The reason is that it is the sum of the water in the ocean plus ice that should be balanced when running ocean-ice models, not the water in any one sub-component. As an extreme example to illustrate the issue, consider an ocean-ice model with zero initial sea ice. As the ocean-ice model spins up, there should be a net accumulation of water in the growing sea ice, and thus a net loss of water from the ocean. The total water contained in the ocean plus ice system is constant, but there is an exchange of water between the subcomponents. This exchange should not be part of the normalization used to balance the hydrological cycle in ocean-ice models. 
+
+
+

trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

@@ -14 +14 @@
 %STEVEN :  is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below
 
-\newpage
-$\ $\newline    % force a new ligne
+%\newpage
+\vspace{2.cm}
+%$\ $\newline    % force a new ligne
 
 Using the representation described in Chap.~\ref{DOM}, several semi-discrete 
@@ -37 +38 @@
 Bottom Boundary Condition), the contribution from the bottom boundary Layer 
 (BBL) parametrisation, and an internal damping (DMP) term. The terms QSR, 
-BBC, BBL and DMP are optional. The external forcings and parameterizations 
+BBC, BBL and DMP are optional. The external forcings and parameterisations 
 require complex inputs and complex calculations (e.g. bulk formulae, estimation 
 of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and 
@@ -54 +55 @@
 
 The different options available to the user are managed by namelist logical or 
-CPP keys. For each equation term ttt, the namelist logicals are \textit{ln\_trattt\_xxx}, 
+CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx}, 
 where \textit{xxx} is a 3 or 4 letter acronym accounting for each optional scheme. 
 The CPP key (when it exists) is \textbf{key\_trattt}. The corresponding code can be 
@@ -62 +63 @@
 equation for output (\key{trdtra} is defined), as described in Chap.~\ref{MISC}.
 
+$\ $\newline    % force a new ligne
 % ================================================================
 % Tracer Advection
@@ -75 +77 @@
 fluxes. Its discrete expression is given by :
 \begin{equation} \label{Eq_tra_adv}
-ADV_\tau =-\frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\left( 
-{\;\delta _i \left[ {e_{2u} {\kern 1pt}e_{3u} {\kern 1pt}\;u\;\tau _u } 
-\right]+\delta _j \left[ {e_{1v} {\kern 1pt}e_{3v} {\kern 1pt}v\;\tau _v } 
-\right]\;} \right)-\frac{1}{\mathop e\nolimits_{3T} }\delta _k \left[ 
-{w\;\tau _w } \right]
-\end{equation}
-where $\tau$ is either T or S. In pure $z$-coordinate (\key{zco} is defined), 
-it reduces to :
+ADV_\tau =-\frac{1}{b_T} \left( 
+\;\delta _i \left[ e_{2u}\,e_{3u} \;  u\; \tau _u  \right]
++\delta _j \left[ e_{1v}\,e_{3v}  \;  v\; \tau _v  \right] \; \right)
+-\frac{1}{e_{3T}} \;\delta _k \left[ w\; \tau _w \right]
+\end{equation}
+where $\tau$ is either T or S, and $b_T= e_{1T}\,e_{2T}\,e_{3T}$ is the volume of $T$-cells. 
+In pure $z$-coordinate (\key{zco} is defined), it reduces to :
 \begin{equation} \label{Eq_tra_adv_zco}
 ADV_\tau =-\frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}}\left( {\;\delta _i 
@@ -89 +90 @@
 e\nolimits_{3T} }\delta _k \left[ {w\;\tau _w } \right]
 \end{equation}
-since the vertical scale factors are functions of $k$ only, and thus $e_{3u} 
-=e_{3v} =e_{3T} $.
-
-The flux form in \eqref{Eq_tra_adv} requires implicitly the use of the continuity equation: 
-$\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
-(using $\nabla \cdot \vect{U}=0)$ or $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$
- in variable volume case ($i.e.$ \key{vvl} defined). Therefore it is of 
-paramount importance to design the discrete analogue of the advection 
-tendency so that it is consistent with the continuity equation in order to 
+since the vertical scale factors are functions of $k$ only, and thus 
+$e_{3u} =e_{3v} =e_{3T} $. The flux form in \eqref{Eq_tra_adv} 
+requires implicitly the use of the continuity equation. Indeed, it is obtained
+by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ 
+which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or 
+$\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant (default option) 
+or variable (\key{vvl} defined) volume case, respectively. 
+Therefore it is of paramount importance to design the discrete analogue of the 
+advection tendency so that it is consistent with the continuity equation in order to 
 enforce the conservation properties of the continuous equations. In other words, 
 by substituting $\tau$ by 1 in (\ref{Eq_tra_adv}) we recover the discrete form of 
@@ -192 +193 @@
 produce a sensible solution. The associated time-stepping is performed using 
 a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in 
-(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
+(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value. The centered second 
+order advection is computed in the \mdl{traadv\_cen2} module. In this module,
+it is also proposed to combine the \textit{cen2} scheme with an upstream scheme
+in specific areas which requires a strong diffusion in order to avoid the generation 
+of false extrema. These areas are the vicinity of large river mouths, some straits 
+with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean 
+temperature is close to the freezing point).
 
 Note that using the cen2 scheme, the overall tracer advection is of second 
 order accuracy since both (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_cen2}) 
-have this order of accuracy.
+have this order of accuracy. Note also that 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -223 +230 @@
 
 A direct consequence of the pseudo-fourth order nature of the scheme is that 
-it is not non-diffusive, i.e. the global variance of a tracer is not 
-preserved using \textit{cen4}. Furthermore, it must be used in conjunction with an 
-explicit diffusion operator to produce a sensible solution. The 
-time-stepping is also performed using a leapfrog scheme in conjunction with 
-an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
+it is not non-diffusive, i.e. the global variance of a tracer is not preserved using 
+\textit{cen4}. Furthermore, it must be used in conjunction with an explicit 
+diffusion operator to produce a sensible solution. The time-stepping is also 
+performed using a leapfrog scheme in conjunction with an Asselin time-filter, 
+so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
 
 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an 
@@ -244 +251 @@
 
 In the Total Variance Dissipation (TVD) formulation, the tracer at velocity 
-points is evaluated using a combination of an upstream and a centred scheme. For 
-example, in the $i$-direction :
+points is evaluated using a combination of an upstream and a centred scheme. 
+For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_tvd}
 \begin{split}
@@ -256 +263 @@
 \end{split}
 \end{equation}
-where $c_u$ is a flux limiter function taking values between 0 and 1. There 
-exist many ways to define $c_u$, each correcponding to a different total 
-variance decreasing scheme. The one chosen in \NEMO is described in 
+where $c_u$ is a flux limiter function taking values between 0 and 1. 
+There exist many ways to define $c_u$, each correcponding to a different 
+total variance decreasing scheme. The one chosen in \NEMO is described in 
 \citet{Zalesak1979}. $c_u$ only departs from $1$ when the advective term 
 produces a local extremum in the tracer field. The resulting scheme is quite 
@@ -264 +271 @@
 This scheme is tested and compared with MUSCL and the MPDATA scheme in 
 \citet{Levy2001}; note that in this paper it is referred to as "FCT" (Flux corrected 
-transport) rather than TVD.
+transport) rather than TVD. The TVD scheme is computed in the \mdl{traadv\_tvd} module.
 
 For stability reasons (see \S\ref{DOM_nxt}), in (\ref{Eq_tra_adv_tvd}) 
@@ -302 +309 @@
 (\np{ln\_traadv\_muscl2}=.true.). Note that the latter choice does not ensure 
 the \textit{positive} character of the scheme. Only the former can be used 
-on both active and passive tracers.
+on both active and passive tracers. The two MUSCL schemes are computed 
+in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -325 +333 @@
 
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 
-error \citep{Sacha2005}. The overall performance of the 
-advection scheme is similar to that reported in \cite{Farrow1995}. 
-It is a relatively good compromise between accuracy and smoothness. It is 
-not a \emph{positive} scheme, meaning that false extrema are permitted, but the 
-amplitude of such are significantly reduced over the centred second order 
-method. Nevertheless it is not recommended that it should be applied to a passive 
-tracer that requires positivity. 
+error \citep{Sacha2005}. The overall performance of the advection 
+scheme is similar to that reported in \cite{Farrow1995}. 
+It is a relatively good compromise between accuracy and smoothness. 
+It is not a \emph{positive} scheme, meaning that false extrema are permitted, 
+but the amplitude of such are significantly reduced over the centred second 
+order method. Nevertheless it is not recommended that it should be applied 
+to a passive tracer that requires positivity. 
 
 The intrinsic diffusion of UBS makes its use risky in the vertical direction 
 where the control of artificial diapycnal fluxes is of paramount importance. 
-Therefore the vertical flux is evaluated using the TVD 
-scheme when \np{ln\_traadv\_ubs}=.true..
+Therefore the vertical flux is evaluated using the TVD scheme when 
+\np{ln\_traadv\_ubs}=.true..
 
 For stability reasons  (see \S\ref{DOM_nxt}), in \eqref{Eq_tra_adv_ubs}, 
@@ -343 +351 @@
 second term (which is the diffusive part of the scheme), is 
 evaluated using the \textit{before} tracer (forward in time). 
-This is discussed by \citet{Webb1998} in the context of the Quick 
-advection scheme. UBS and QUICK 
-schemes only differ by one coefficient. Replacing 1/6 with 1/8 in 
-\eqref{Eq_tra_adv_ubs} leads to the QUICK advection scheme 
-\citep{Webb1998}. This option is not available through a namelist 
-parameter, since the 1/6 coefficient is hard coded. Nevertheless 
-it is quite easy to make the substitution in the \mdl{traadv\_ubs} module 
-and obtain a QUICK scheme.
+This choice is discussed by \citet{Webb1998} in the context of the 
+QUICK advection scheme. UBS and QUICK schemes only differ 
+by one coefficient. Replacing 1/6 with 1/8 in \eqref{Eq_tra_adv_ubs} 
+leads to the QUICK advection scheme \citep{Webb1998}. 
+This option is not available through a namelist parameter, since the 
+1/6 coefficient is hard coded. Nevertheless it is quite easy to make the 
+substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
 
 Note that :
 
-(1): When a high vertical resolution $O(1m)$ is used, the model stability can 
+(1) When a high vertical resolution $O(1m)$ is used, the model stability can 
 be controlled by vertical advection (not vertical diffusion which is usually 
 solved using an implicit scheme). Computer time can be saved by using a 
-time-splitting technique on vertical advection. This case has been 
-implemented and validated in ORCA05 with 301 levels. It is not available in the 
-current reference version. 
-
-(2) : In a forthcoming release four options will be available for the vertical 
+time-splitting technique on vertical advection. Such a technique has been 
+implemented and validated in ORCA05 with 301 levels. It is not available 
+in the current reference version. 
+
+(2) In a forthcoming release four options will be available for the vertical 
 component used in the UBS scheme. $\tau _w^{ubs}$ will be evaluated 
-using either \textit{(a)} a centred $2^{nd}$ order scheme , or  \textit{(b)} 
+using either \textit{(a)} a centred $2^{nd}$ order scheme, or  \textit{(b)} 
 a TVD scheme, or  \textit{(c)} an interpolation based on conservative 
 parabolic splines following the \citet{Sacha2005} implementation of UBS 
@@ -369 +376 @@
 similar to an eighth-order accurate conventional scheme.
 
-following \citet{Sacha2005} implementation of UBS in ROMS, or  \textit{(d)} 
-an UBS. The $3^{rd}$ case has dispersion properties similar to an 
-eight-order accurate conventional scheme.
-
-(3) : It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
-\begin{equation} \label{Eq_tra_adv_ubs2}
-\tau _u^{ubs} = \left\{  \begin{aligned}
-   & \tau _u^{cen4} + \frac{1}{12} \tau"_i      & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
-   & \tau _u^{cen4} - \frac{1}{12} \tau"_{i+1}  & \quad \text{if }\ u_{i+1/2}       <       0
-                   \end{aligned}    \right.
+(3) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
+\begin{equation} \label{Eq_traadv_ubs2}
+\tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{  
+   \begin{aligned}
+   & + \tau"_i       & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\
+   &  - \tau"_{i+1}     & \quad \text{if }\ u_{i+1/2}       <       0
+   \end{aligned}    \right.
 \end{equation}
 or equivalently 
-\begin{equation} \label{Eq_tra_adv_ubs2b}
+\begin{equation} \label{Eq_traadv_ubs2b}
 u_{i+1/2} \ \tau _u^{ubs} 
 =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta _i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2}
 - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
 \end{equation}
-\eqref{Eq_tra_adv_ubs2} has several advantages. Firstly, it clearly reveals 
+
+\eqref{Eq_traadv_ubs2} has several advantages. Firstly, it clearly reveals 
 that the UBS scheme is based on the fourth order scheme to which an 
 upstream-biased diffusion term is added. Secondly, this emphasises that the 
@@ -394 +399 @@
 coefficient which is simply proportional to the velocity:
  $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v2.3 still uses 
- \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_tra_adv_ubs2}. This should be 
+ \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}. This should be 
  changed in forthcoming release.
  %%%
@@ -411 +416 @@
 is the third order Godunov scheme. It is associated with the ULTIMATE QUICKEST 
 limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray 
-(MERCATOR-ocean). 
-The resulting scheme is quite expensive but \emph{positive}. It can be used on both active and passive tracers. Nevertheless, the intrinsic diffusion of QCK makes its use 
-risky in the vertical direction where the control of artificial diapycnal fluxes is of 
-paramount importance. Therefore the vertical flux is evaluated using the CEN2 
-scheme. This no more ensure the positivity of the scheme. The use of TVD in the 
-vertical direction as for the UBS case should be implemented to maintain the property.
+(MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
+The resulting scheme is quite expensive but \emph{positive}. 
+It can be used on both active and passive tracers. 
+Nevertheless, the intrinsic diffusion of QCK makes its use risky in the vertical 
+direction where the control of artificial diapycnal fluxes is of paramount importance. 
+Therefore the vertical flux is evaluated using the CEN2 scheme. 
+This no more ensure the positivity of the scheme. The use of TVD in the vertical 
+direction as for the UBS case should be implemented to maintain the property.
 
 
@@ -430 +437 @@
 with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented 
 in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference 
-version 2.3.
+version 3.0.
 
 % ================================================================
@@ -447 +454 @@
 coefficients (either constant or variable in space and time) as well as the 
 computation of the slope along which the operators act, are performed in the 
-\mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. The lateral diffusion of tracers is evaluated using a forward scheme, 
+\mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. 
+The lateral diffusion of tracers is evaluated using a forward scheme, 
 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time, 
 except for the pure vertical component that appears when a rotation tensor 
@@ -456 +464 @@
 %        Iso-level laplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Iso-level laplacian operator (\textit{traldf\_lap} - \np{ln\_traldf\_lap})]
-         {Iso-level laplacian operator (\mdl{traldf\_lap} - \np{ln\_traldf\_lap}=.true.) }
+\subsection   [Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap})]
+         {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=.true.) }
 \label{TRA_ldf_lap}
 
-A laplacian diffusion operator (i.e. a harmonic operator) acting along the model 
+A laplacian diffusion operator ($i.e.$ a harmonic operator) acting along the model 
 surfaces is given by: 
 \begin{equation} \label{Eq_tra_ldf_lap}
-\begin{split}
-D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\;  e_{3T} } &\left[ {\quad \delta _i 
-\left[ {A_u^{lT} \left( {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} 
-\left[ T \right]} \right)} \right]} \right.
-\\
-&\ \left. {+\; \delta _j \left[ 
-{A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ T 
-\right]} \right)} \right]\quad } \right]
-\end{split}
-\end{equation}
-
-This lateral operator is a \emph{horizontal} one ($i.e.$ acting along 
-geopotential surfaces) in the $z$-coordinate with or without partial step, 
-but is simply an iso-level operator in the $s$-coordinate. 
+D_T^{lT} =\frac{1}{b_T} \left( \;
+   \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right] 
++ \delta _{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right]  \;\right)
+\end{equation}
+where  $b_T= e_{1T}\,e_{2T}\,e_{3T}$  is the volume of $T$-cells. 
+It can be found in the \mdl{traadv\_lap} module.
+
+This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal} 
+operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate with 
+or without partial step, but is simply an iso-level operator in the $s$-coordinate. 
 It is thus used when, in addition to \np{ln\_traldf\_lap}=.true., we have 
-\np{ln\_traldf\_level}=.true., or both \np{ln\_traldf\_hor}=.true. and 
-\np{ln\_zco}=.false.. In both cases, it significantly contributes to 
-diapycnal mixing. It is therefore not recommended.
+\np{ln\_traldf\_level}=.true., or \np{ln\_traldf\_hor}=\np{ln\_zco}=.true.. 
+In both cases, it significantly contributes to diapycnal mixing. 
+It is therefore not recommended.
 
 Note that 
-(1) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}=e_{3v}=e_{3T}$, so 
-that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}). 
-(2) In partial step $z$-coordinate (\np{ln\_zps}=.true.), tracers in horizontally 
+(a) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}=e_{3v}=e_{3T}$, 
+so that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}) ; 
+(b) In partial step $z$-coordinate (\np{ln\_zps}=.true.), tracers in horizontally 
 adjacent cells are located at different depths in the vicinity of the bottom. 
 In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level 
@@ -494 +498 @@
 %        Rotated laplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Rotated laplacian operator (\textit{traldf\_iso} - \np{ln\_traldf\_lap})]
-         {Rotated laplacian operator (\mdl{traldf\_iso} - \np{ln\_traldf\_lap}=.true.)}
+\subsection   [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})]
+         {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=.true.)}
 \label{TRA_ldf_iso}
 
@@ -501 +505 @@
 (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and 
 $s$-coordinates:
-
 \begin{equation} \label{Eq_tra_ldf_iso}
 \begin{split}
- D_T^{lT} =& \frac{1}{e_{1T}\,e_{2T}\,e_{3T} }
- \\
-& \left\{ {\delta _i \left[ {A_u^{lT}  \left( 
-    {\frac{e_{2u} \; e_{3u} }{e_{1u} } \,\delta _{i+1/2}[T]
-   -e_{2u} \; r_{1u} \,\overline{\overline {\delta _{k+1/2}[T]}}^{\,i+1/2,k}}
- \right)} \right]} \right. 
-\\ 
-& +\delta 
-_j \left[ {A_v^{lT} \left( {\frac{e_{1v}\,e_{3v} }{e_{2v} 
-}\,\delta _{j+1/2} \left[ T \right]-e_{1v}\,r_{2v} 
-\,\overline{\overline {\delta _{k+1/2} \left[ T \right]}} ^{\,j+1/2,k}} 
-\right)} \right] 
-\\ 
-& +\delta 
-_k \left[ {A_w^{lT} \left( 
--e_{2w}\,r_{1w} \,\overline{\overline {\delta _{i+1/2} \left[ T \right]}} ^{\,i,k+1/2}
-\right.} \right. 
-\\ 
+ D_T^{lT} = \frac{1}{b_T}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left( 
+     \frac{e_{2u}\;e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T]
+   - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k}
+                                                     \right)   \right]   \right.    \\ 
+&             +\delta_j \left[ A_v^{lT} \left( 
+          \frac{e_{1v}\,e_{3v}}{e_{2v}}  \,\delta_{j+1/2} [T] 
+        - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} 
+                                                    \right)   \right]                 \\ 
+& +\delta_k \left[ A_w^{lT} \left( 
+       -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2}
+                                                    \right.   \right.                 \\ 
 & \qquad \qquad \quad 
--e_{1w}\,r_{2w} \,\overline{\overline {\delta _{j+1/2} \left[ T \right]}} ^{\,j,k+1/2}
-\\
-& \left. {\left. { 
- \quad \quad \quad \left. {{\kern 
-1pt}+\frac{e_{1w}\,e_{2w} }{e_{3w} }\,\left( {r_{1w} ^2+r_{2w} ^2} 
-\right)\,\delta _{k+1/2} \left[ T \right]} \right)} \right]\;\;\;} \right\} 
+        - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2}     \\
+& \left. {\left. {   \qquad \qquad \ \ \ \left. {
+        +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right)
+           \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 
  \end{split}
  \end{equation}
-where $r_1$ and $r_2$ are the slopes between the surface of computation 
+where $b_T= e_{1T}\,e_{2T}\,e_{3T}$  is the volume of $T$-cells, 
+$r_1$ and $r_2$ are the slopes between the surface of computation 
 ($z$- or $s$-surfaces) and the surface along which the diffusion operator 
 acts ($i.e.$ horizontal or iso-neutral surfaces).  It is thus used when, 
-in addition to \np{ln\_traldf\_lap}=.true., we have \np{ln\_traldf\_iso}=.true., 
+in addition to \np{ln\_traldf\_lap}= .true., we have \np{ln\_traldf\_iso}=.true., 
 or both \np{ln\_traldf\_hor}=.true. and \np{ln\_zco}=.true.. The way these 
 slopes are evaluated is given in \S\ref{LDF_slp}. At the surface, bottom 
@@ -544 +539 @@
 be solved using the same implicit time scheme as that used in the vertical 
 physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term 
-is not computed in the \mdl{traldf} module, but in the \mdl{trazdf} module 
+is not computed in the \mdl{traldf\_iso} module, but in the \mdl{trazdf} module 
 where, if iso-neutral mixing is used, the vertical mixing coefficient is simply 
 increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$. 
@@ -552 +547 @@
 (see \S\ref{LDF}) allows the model to run safely without any additional 
 background horizontal diffusion \citep{Guily2001}. An alternative scheme 
-\citep{Griffies1998} which preserves both tracer and its variance is currently 
-been tested in \NEMO. 
+developed by \cite{Griffies1998} which preserves both tracer and its variance 
+is currently been tested in \NEMO. It should be available in a forthcoming
+release.
 
 Note that in the partial step $z$-coordinate (\np{ln\_zps}=.true.), the horizontal 
@@ -562 +558 @@
 %        Iso-level bilaplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Iso-level bilaplacian operator (\textit{traldf\_bilap} - \np{ln\_traldf\_bilap})]
-         {Iso-level bilaplacian operator (\mdl{traldf\_bilap} - \np{ln\_traldf\_bilap}=.true.)}
+\subsection   [Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap})]
+         {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=.true.)}
 \label{TRA_ldf_bilap}
 
@@ -569 +565 @@
 applying (\ref{Eq_tra_ldf_lap}) twice. It requires an additional assumption 
 on boundary conditions: the first and third derivative terms normal to the 
-coast are set to zero.
-
-It is used when, in addition to \np{ln\_traldf\_bilap}=.true., we have 
-\np{ln\_traldf\_level}=.true., or both \np{ln\_traldf\_hor}=.true. and 
+coast are set to zero. It is used when, in addition to \np{ln\_traldf\_bilap}=.true., 
+we have \np{ln\_traldf\_level}=.true., or both \np{ln\_traldf\_hor}=.true. and 
 \np{ln\_zco}=.false.. In both cases, it can contribute diapycnal mixing, 
 although less than in the laplacian case. It is therefore not recommended.
 
 Note that in the code, the bilaplacian routine does not call the laplacian 
-routine twice but is rather a separate routine. This is due to the fact that we 
-introduce the eddy diffusivity coefficient, A, in the operator as: $\nabla 
-\cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$, instead of 
-$-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ where 
-$a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations 
-ensure the total variance decrease, but the former requires a larger number 
-of code-lines. It will be corrected in a forthcoming release.
+routine twice but is rather a separate routine that can be found in the
+\mdl{traldf\_bilap} module. This is due to the fact that we introduce the 
+eddy diffusivity coefficient, A, in the operator as: 
+$\nabla \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$, 
+instead of 
+$-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ 
+where $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations 
+ensure the total variance decrease, but the former requires a larger 
+number of code-lines. It will be corrected in a forthcoming release.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Rotated bilaplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Rotated bilaplacian operator (\textit{traldf\_bilapg} - \np{ln\_traldf\_bilap})]
-         {Rotated bilaplacian operator (\mdl{traldf\_bilapg} - \np{ln\_traldf\_bilap}=.true.)}
+\subsection   [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})]
+         {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=.true.)}
 \label{TRA_ldf_bilapg}
 
@@ -595 +591 @@
 applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption 
 on boundary conditions: first and third derivative terms normal to the 
-coast, the bottom and the surface are set to zero.
-
-It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have 
-\np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T. 
+coast, the bottom and the surface are set to zero. It can be found in the
+\mdl{traldf\_bilapg}.
+
+It is used when, in addition to \np{ln\_traldf\_bilap}=.true., we have 
+\np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=.true. and \np{ln\_zco}=.true.. 
 Nevertheless, this rotated bilaplacian operator has never been seriously 
 tested. No warranties that it is neither free of bugs or correctly formulated. 
@@ -619 +616 @@
 The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the 
 following semi-discrete space form:
-(\ref{Eq_PE_zdf}) takes the following semi-discrete space form:
 \begin{equation} \label{Eq_tra_zdf}
 \begin{split}
-D^{vT}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ 
-\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T]   \right] 
+D^{vT}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right] 
 \\
-D^{vS}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ 
-\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S]   \right] 
+D^{vS}_T &= \frac{1}{e_{3T}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S] \;\right] 
 \end{split}
 \end{equation}
-
 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity 
-coefficients on Temperature and Salinity, respectively. Generally, 
+coefficients on temperature and salinity, respectively. Generally, 
 $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is 
-parameterised (\key{zdfddm} is defined). The way these coefficients 
+parameterised ($i.e.$ \key{zdfddm} is defined). The way these coefficients 
 are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when 
 iso-neutral mixing is used, both mixing coefficients are increased 
@@ -640 +633 @@
 
 At the surface and bottom boundaries, the turbulent fluxes of 
-momentum, heat and salt must be specified. At the surface they 
-are prescribed from the surface forcing (see \S\ref{TRA_sbc}), 
+heat and salt must be specified. At the surface they are prescribed 
+from the surface forcing and added in a dedicated routine (see \S\ref{TRA_sbc}), 
 whilst at the bottom they are set to zero for heat and salt unless 
 a geothermal flux forcing is prescribed as a bottom boundary 
-condition (\S\ref{TRA_bbc}). 
+condition (see \S\ref{TRA_bbc}). 
 
 The large eddy coefficient found in the mixed layer together with high 
@@ -712 +705 @@
  \end{aligned}
 \end{equation} 
-
 where EMP is the freshwater budget (evaporation minus precipitation 
 minus river runoff) which forces the ocean volume, $Q_{ns}$ is the 
@@ -722 +714 @@
 (except for the effect of the Asselin filter).
 
-%AMT note: the ice-ocean flux had been forgotten in the first release of the key_vvl option, has this been corrected in the code? 
+%AMT note: the ice-ocean flux had been forgotten in the first release of the key_vvl option, has this been corrected in the code?     ===> gm :  NO to be added at NOCS 
 
 In the second case (linear free surface), the vertical scale factors are 
@@ -729 +721 @@
 temperature of precipitation and runoffs, $F_{wf}^e +F_{wf}^i =0$ 
 for temperature. The resulting forcing term for temperature is: 
-
 \begin{equation} \label{Eq_tra_forcing_q}
 F^T=\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} }
@@ -779 +770 @@
 \end{equation}
 
-where $I$ is the downward irradiance. The additional term in \eqref{Eq_PE_qsr} is discretized as follows:
+where $I$ is the downward irradiance. The additional term in \eqref{Eq_PE_qsr} 
+is discretized as follows:
 \begin{equation} \label{Eq_tra_qsr}
 \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3T}} \delta_k \left[ I_w \right]
@@ -796 +788 @@
 \gmcomment : Jerlov reference to be added
 %
-classification: $\xi_1 = 0.35m$, $\xi_2 = 0.23m$ and $R = 0.58$ 
+classification: $\xi_1 = 0.35~m$, $\xi_2 = 23~m$ and $R = 0.58$ 
 (corresponding to \np{xsi1}, \np{xsi2} and \np{rabs} namelist parameters, 
 respectively). $I$ is masked (no flux through the ocean bottom), 
@@ -837 +829 @@
 \begin{figure}[!t] \label{Fig_geothermal}  \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_geoth.pdf}
-\caption{Geothermal Heat flux (in $mW.m^{-2}$) as inferred from the age 
-of the sea floor and the formulae of \citet{Stein1992}.}
+\caption{Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OSD08}.
+It is inferred from the age of the sea floor and the formulae of \citet{Stein1992}.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -845 +837 @@
 the ocean bottom, $i.e.$ a no flux boundary condition is applied on active 
 tracers at the bottom. This is the default option in \NEMO, and it is 
-implemented using the masking technique. Hoever, there is a 
+implemented using the masking technique. However, there is a 
 non-zero heat flux across the seafloor that is associated with solid 
 earth cooling. This flux is weak compared to surface fluxes (a mean 
 global value of $\sim0.1\;W/m^2$ \citep{Stein1992}), but it is 
-systematically positive and acts on the densest water masses. Taking 
-this flux into account in a global ocean model increases
-the deepest overturning cell (i.e. the one associated with the Antarctic 
-Bottom Water) by a few Sverdrups. 
+systematically positive and acts on the densest water masses. 
+Taking this flux into account in a global ocean model increases
+the deepest overturning cell ($i.e.$ the one associated with the Antarctic 
+Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OSD08}. 
 
 The presence or not of geothermal heating is controlled by the namelist 
@@ -886 +878 @@
 a sill \citep{Willebrand2001}, and the thickness of the plume is not resolved. 
 
-The idea of the bottom boundary layer (BBL) parameterization first introduced by 
+The idea of the bottom boundary layer (BBL) parameterisation first introduced by 
 \citet{BeckDos1998} is to allow a direct communication between 
 two adjacent bottom cells at different levels, whenever the densest water is 
@@ -933 +925 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_bbl}  \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
+\includegraphics[width=0.8\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
 \caption{Advective Bottom Boundary Layer.}
 \end{center}   \end{figure}
@@ -1047 +1039 @@
 \end{split}
 \end{equation} 
-
 where $\text{RHS}_T$ is the right hand side of the temperature equation, 
 the subscript $f$ denotes filtered values and $\gamma$ is the Asselin 
@@ -1093 +1084 @@
 the practical salinity (another \NEMO variable) and the pressure (assuming no 
 pressure variation along geopotential surfaces, i.e. the pressure in decibars is 
-approximated by the depth in meters). Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state have exactly the same except that 
+approximated by the depth in meters). Both the \citet{UNESCO1983} and 
+\citet{JackMcD1995} equations of state have exactly the same except that 
 the values of the various coefficients have been adjusted by \citet{JackMcD1995} 
 in order to directly use the \textit{potential} temperature instead of the 
@@ -1194 +1186 @@
 \begin{equation} \label{Eq_tra_eos_fzp}
    \begin{split}
-T_f (S,p) &= \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
+T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} 
                        -  2.154996 \;10^{-4} \,S  \right) \ S    \\
-               & - 7.53\,10^{-3}\,p 
+               - 7.53\,10^{-3} \ \ p 
    \end{split}
 \end{equation}

trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex

@@ -7 +7 @@
 
 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN.
-\gmcomment{Steven remark : problem here with turbulent vs turbulence.  I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation!
-Gurvan :  I kept "turbulent closure"...}
-\gmcomment{Steven bis : parameterization is the american spelling, parameterisation is the british}
+\gmcomment{Steven remark (not taken into account : problem here with turbulent vs turbulence.  I've changed "turbulent closure" to "turbulence closure", "turbulent mixing length" to "turbulence mixing length", but I've left "turbulent kinetic energy" alone - though I think it is an historical abberation!
+Gurvan :  I kept "turbulent closure etc "...}
 
 
@@ -25 +24 @@
 flux forcing is prescribed as a bottom boundary condition ($i.e.$ \key{trabbl} 
 defined, see \S\ref{TRA_bbc}), and specified through a bottom friction 
-parameterization for momentum (see \S\ref{ZDF_bfr}).
+parameterisation for momentum (see \S\ref{ZDF_bfr}).
 
 In this section we briefly discuss the various choices offered to compute 
@@ -84 +83 @@
 large scale ocean structures. The hypothesis of a mixing mainly maintained by the 
 growth of Kelvin-Helmholtz like instabilities leads to a dependency between the 
-vertical turbulence eddy coefficients and the local Richardson number ($i.e.$ the 
+vertical eddy coefficients and the local Richardson number ($i.e.$ the 
 ratio of stratification to vertical shear). Following \citet{PacPhil1981}, the following 
 formulation has been implemented:
@@ -114 +113 @@
 The vertical eddy viscosity and diffusivity coefficients are computed from a TKE 
 turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent 
-kinetic energy, and a closure assumption for the turbulence length scales. This 
+kinetic energy, and a closure assumption for the turbulent length scales. This 
 turbulent closure model has been developed by \citet{Bougeault1989} in the 
 atmospheric case, adapted by \citet{Gaspar1990} for the oceanic case, and 
@@ -156 +155 @@
 The choice of $P_{rt} $ is controlled by the \np{npdl} namelist parameter.
 
-For computational efficiency, the original formulation of the turbulence length 
+For computational efficiency, the original formulation of the turbulent length 
 scales proposed by \citet{Gaspar1990} has been simplified. Four formulations 
 are proposed, the choice of which is controlled by the \np{nmxl} namelist 
@@ -186 +185 @@
 constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, 
 the upward and downward length scales, and evaluate the dissipation and 
-mixing turbulence length scales as (and note that here we use numerical 
+mixing turbulent length scales as (and note that here we use numerical 
 indexing):
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -208 +207 @@
 In the \np{nmxl}=2 case, the dissipation and mixing length scales take the same 
 value: $ l_k=  l_\epsilon = \min \left(\ l_{up} \;,\;  l_{dwn}\ \right)$, while in the 
-\np{nmxl}=2 case, the dissipation and mixing turbulence length scales are give 
+\np{nmxl}=2 case, the dissipation and mixing length scales are give 
 as in \citet{Gaspar1990}:
 \begin{equation} \label{Eq_tke_mxl_gaspar}
@@ -243 +242 @@
 %--------------------------------------------------------------------------------------------------------------
 
-The KKP scheme has been implemented by J. Chanut ...
+The K-Profile Parametrization (KKP) developed by \cite{Large_al_RG94} has been 
+implemented in \NEMO by J. Chanut (PhD reference to be added here!).
 
 \colorbox{yellow}{Add a description of KPP here.}
@@ -262 +262 @@
 quickly re-establish the static stability of the water column. These 
 processes have been removed from the model via the hydrostatic 
-assumption so they must be parameterized. Three parameterizations 
+assumption so they must be parameterized. Three parameterisations 
 are available to deal with convective processes: a non-penetrative 
 convective adjustment or an enhanced vertical diffusion, or/and the 
@@ -354 +354 @@
 %--------------------------------------------------------------------------------------------------------------
 
-The enhanced vertical diffusion parameterization is used when \np{ln\_zdfevd}=true. 
+The enhanced vertical diffusion parameterisation is used when \np{ln\_zdfevd}=true. 
 In this case, the vertical eddy mixing coefficients are assigned very large values 
 (a typical value is $10\;m^2s^{-1})$ in regions where the stratification is unstable 
@@ -364 +364 @@
 if \np{n\_evdm}=1, the four neighbouring $A_u^{vm} \;\mbox{and}\;A_v^{vm}$ 
 values also, are set equal to the namelist parameter \np{avevd}. A typical value 
-for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterization of 
+for $avevd$ is between 1 and $100~m^2.s^{-1}$. This parameterisation of 
 convective processes is less time consuming than the convective adjustment 
 algorithm presented above when mixing both tracers and momentum in the 
@@ -384 +384 @@
 $A_u^{vm} {and}\;A_v^{vm}$ (up to $1\;m^2s^{-1})$. These large values 
 restore the static stability of the water column in a way similar to that of the 
-enhanced vertical diffusion parameterization (\S\ref{ZDF_evd}). However, 
+enhanced vertical diffusion parameterisation (\S\ref{ZDF_evd}). However, 
 in the vicinity of the sea surface (first ocean layer), the eddy coefficients 
-computed by the turbulence scheme do not usually exceed $10^{-2}m.s^{-1}$, 
+computed by the turbulent closure scheme do not usually exceed $10^{-2}m.s^{-1}$, 
 because the mixing length scale is bounded by the distance to the sea surface 
 (see \S\ref{ZDF_tke}). It can thus be useful to combine the enhanced vertical 
@@ -412 +412 @@
 to diffusive convection. Double-diffusive phenomena contribute to diapycnal 
 mixing in extensive regions of the ocean.  \citet{Merryfield1999} include a 
-parameterization of such phenomena in a global ocean model and show that 
+parameterisation of such phenomena in a global ocean model and show that 
 it leads to relatively minor changes in circulation but exerts significant regional 
 influences on temperature and salinity. 
@@ -422 +422 @@
 \end{align*}
 where subscript $f$ represents mixing by salt fingering, $d$ by diffusive convection, 
-and $o$ by processes other than double diffusion. The rates of double-diffusive mixing depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 
+and $o$ by processes other than double diffusion. The rates of double-diffusive mixing 
+depend on the buoyancy ratio $R_\rho = \alpha \partial_z T / \beta \partial_z S$, 
 where $\alpha$ and $\beta$ are coefficients of thermal expansion and saline 
 contraction (see \S\ref{TRA_eos}). To represent mixing of $S$ and $T$ by salt 
@@ -443 +444 @@
 $A^{\ast v} = 10^{-4}~m^2.s^{-1}$ ; (b) diapycnal diffusivities $A_d^{vT}$ and 
 $A_d^{vS}$ for temperature and salt in regions of diffusive convection. Heavy 
-curves denote the Federov parameterization and thin curves the Kelley 
-parameterization. The latter is not implemented in \NEMO. }
+curves denote the Federov parameterisation and thin curves the Kelley 
+parameterisation. The latter is not implemented in \NEMO. }
 \end{center}    \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -453 +454 @@
 we adopt $R_c = 1.6$, $n = 6$, and $A^{\ast v} = 10^{-4}~m^2.s^{-1}$.
 
-To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested by Federov (1988) is used: 
+To represent mixing of S and T by diffusive layering,  the diapycnal diffusivities suggested 
+by Federov (1988) is used: 
 \begin{align}  \label{Eq_zdfddm_d}
 A_d^{vT} &=    \begin{cases}
@@ -525 +527 @@
 \label{ZDF_bfr_linear}
 
-The linear bottom friction parameterization assumes that the bottom friction 
+The linear bottom friction parameterisation assumes that the bottom friction 
 is proportional to the interior velocity (i.e. the velocity of the last model level):
 \begin{equation} \label{Eq_zdfbfr_linear}
@@ -571 +573 @@
 \label{ZDF_bfr_nonlinear}
 
-The non-linear bottom friction parameterization assumes that the bottom 
+The non-linear bottom friction parameterisation assumes that the bottom 
 friction is quadratic: 
 \begin{equation} \label{Eq_zdfbfr_nonlinear}

trunk/DOC/TexFiles/Chapters/Introduction.tex

@@ -15 +15 @@
 its interactions with the other components of the earth climate system (atmosphere, 
 sea-ice, biogeochemical tracers, ...) over a wide range of space and time scales. 
-This documentation provides information about the physics represented by the ocean component of \NEMO and the rationale for the choice of numerical schemes and 
+This documentation provides information about the physics represented by the ocean 
+component of \NEMO and the rationale for the choice of numerical schemes and 
 the model design. More specific information about running the model on different 
 computers, or how to set up a configuration, are found on the \NEMO web site 
@@ -92 +93 @@
 around the code, the module names follow a three-letter rule. For example, \mdl{tradmp} 
 is a module related to the TRAcers equation, computing the DaMPing. The complete list 
-of module names is presented in \colorbox{yellow}{annex}. Furthermore, modules are  
+of module names is presented in Appendix~\ref{Apdx_D}. Furthermore, modules are  
 organized in a few directories that correspond to their category, as indicated by the first 
 three letters of their name. 
@@ -115 +116 @@
 \end{table}
 
-In the current release (v2.3), LBC directory (see Chap.~\ref{LBC}) does not yet exist. 
-When created LBC will contain the OBC directory (Open Boundary Condition), and the
-\mdl{lbclnk}, \mdl{mppini} and \mdl{lib\_mpp} modules. 
+In the current release (v3.0), the LBC directory does not yet exist. 
+When created LBC will contain the OBC directory (Open Boundary Condition), 
+and the \mdl{lbclnk}, \mdl{mppini} and \mdl{lib\_mpp} modules. 
 
  \vspace{1cm}   Nota Bene : \vspace{0.25cm}
@@ -141 +142 @@
 (9) online diagnostics : tracers trend in the mixed layer and vorticity balance; \\
 (10) rewriting of the I/O management; \\
-(11) OASIS 3 and 4 couplers interfacing with atmospheric global circulation models. 
+(11) OASIS 3 and 4 couplers interfacing with atmospheric global circulation models. \\
 (12) surface module (SBC) that simplify the way the ocean is forced and include two
-bulk formulea (CLIO and CORE)
+bulk formulea (CLIO and CORE)\\
 (13) introduction of LIM 3, the new Louvain-la-Neuve sea-ice model (C-grid rheology and
-new thermodynamics including bulk ice salinity)
+new thermodynamics including bulk ice salinity) \citep{Vancoppenolle_al_OM08}
 
 In addition, several minor modifications in the coding have been introduced with the constant concern of improving performance on both scalar and vector computers.