Changeset 817 for trunk/DOC/BETA/Chapters/Chap_DOM.tex

Timestamp:

2008-02-09T15:13:48+01:00 (16 years ago)

Author:

gm

Message:

trunk - update including Steven correction of the first 5 chapters (until DYN) and activation of Appendix A & B

File:

• trunk/DOC/BETA/Chapters/Chap_DOM.tex (modified) (34 diffs)

trunk/DOC/BETA/Chapters/Chap_DOM.tex

@@ -8 +8 @@
 
 % Missing things:
-%  - istate: description of the initial state
+%  - istate: description of the initial state   ==> this has to be put elsewhere..
+%                  perhaps in MISC ?  By the way the initialisation of T S and dynamics 
+%                  should be put outside of DOM routine (better with TRC staff and off-line
+%                  tracers)
 %  - daymod: definition of the time domain (nit000, nitend andd the calendar)
 %  -geo2ocean:  how to switch from geographic to mesh coordinate
 %  - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled
 
-
-
-
-
-Having defined the continuous equations in Chap.~\ref{PE}, we need to choose a discretization on a grid, and numerical algorithms. In the present chapter, we provide a general description of the staggered grid used in OPA, and other information relevant to the main directory routines (time stepping, main program) as well as the DOM (DOMain) directory. 
+\gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, would help  ==> to be added}
+%%%%
+
+
+\newpage
+$\ $\newline    % force a new ligne
+
+
+Having defined the continuous equations in Chap.~\ref{PE}, we need to choose a 
+discretization on a grid, and numerical algorithms. In the present chapter, we 
+provide a general description of the staggered grid used in \NEMO, and other 
+information relevant to the main directory routines (time stepping, main program) 
+as well as the DOM (DOMain) directory. 
 
 % ================================================================
@@ -34 +45 @@
 \begin{figure}[!tb] \label{Fig_cell}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./Figures/Fig_cell.pdf}
-\caption{Arrangement of variables. $T$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) indicates vector points, and $f$ indicates vorticity points where both relative and planetary vorticities are defined}
+\caption{Arrangement of variables. $T$ indicates scalar points where temperature, 
+salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) 
+indicates vector points, and $f$ indicates vorticity points where both relative and 
+planetary vorticities are defined}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. Special attention has been given to the homogeneity of the solution in the three space directions. The arrangement of variables is the 
-same in all directions. It consists in cells centred on scalar points ($T$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.
-
-The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. Each scale factor is defined as the local analytical value 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 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 scale factors must be evaluated by centred finite difference approximation, not from their 
-analytical expression. This preserves the symmetry of the discrete set of equations and therefore allows satisfying many of the continuous properties (see { \colorbox{yellow}{Annexe C}). A similar, related remark can be made about the domain size: when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). 
+The numerical techniques used to solve the Primitive Equations in this model are 
+based on the traditional, centred second-order finite difference approximation. 
+Special attention has been given to the homogeneity of the solution in the three 
+space directions. The arrangement of variables is the same in all directions. 
+It consists of cells centred on scalar points ($T$, $S$, $p$, $\rho$) with vector 
+points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). 
+This is the generalisation to three dimensions of the well-known ``C'' grid in 
+Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and 
+planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge 
+and the barotropic stream function $\psi$ is defined at horizontal points overlying 
+the $\zeta$ and $f$-points.
+
+The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined 
+by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. 
+The grid-points are located at integer or integer and a half value of $(i,j,k)$ as 
+indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$, 
+$f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale 
+factors are defined. Each scale factor is defined as the local analytical value 
+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 
+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 
+scale factors must be evaluated by centred finite difference approximation, not 
+from their analytical expression. This preserves the symmetry of the discrete set 
+of equations and therefore satisfies many of the continuous properties (see 
+Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain 
+size: when needed, an area, volume, or the total ocean depth must be evaluated 
+as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). 
 
 \begin{table}[!tb] \label{Tab_cell}
@@ -56 +94 @@
 fw & $i+1/2$   & $j+1/2$   & $k+1/2$   \\ \hline
 \end{tabular}
-\caption{Location of grid-points as a function of integer or integer and a half value of the column, line or level. This indexation is only used for the writing of semi-discrete equation. In the code, the indexation use integer value only and has a reverse direction in the vertical (see \S\ref{DOM_Num_Index})}
+\caption{Location of grid-points as a function of integer or integer and a half value 
+of the column, line or level. This indexing is only used for the writing of the semi-
+discrete equation. In the code, the indexing uses integer values only and has a 
+reverse direction in the vertical (see \S\ref{DOM_Num_Index})}
 \end{center}
 \end{table}
@@ -66 +107 @@
 \label{DOM_operators}
 
-Given the values of a variable $q$ at adjacent points, the derivation and averaging operators at the midpoint between them are:
+Given the values of a variable $q$ at adjacent points, the differencing and 
+averaging operators at the midpoint between them are:
 \begin{subequations} \label{Eq_di_mi}
 \begin{align}
- \delta _i [q]  &=  \  \    q(i+1/2)  - q(i-1/2)      \\
- \overline q^i &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2
+ \delta _i [q]       &=  \  \    q(i+1/2)  - q(i-1/2)    \\
+ \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2
 \end{align}
 \end{subequations}
 
-Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$. Following 
-\eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a variable $q$ defined at $T$-point has its three components defined at $(u,v,w)$ while its Laplacien is defined at $T$-point. These operators have the following discrete forms in the curvilinear $s$-coordinate system:
+Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and 
+$k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a 
+variable $q$ defined at a $T$-point has its three components defined at $u$-, $v$- 
+and $w$-points while its Laplacien is defined at $T$-point. These operators have 
+the following discrete forms in the curvilinear $s$-coordinate system:
 \begin{equation} \label{Eq_DOM_grad}
 \nabla q\equiv    \frac{1}{e_{1u} }\delta _{i+1/2} \left[ q \right]\;\,{\rm {\bf i}}
@@ -91 +136 @@
 \end{multline}
 
-Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl components defined at $(vw,uw,f)$ and its divergence defined at $T$-points:
+Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl 
+components defined at $vw$-, $uw$, and $f$-points, and its divergence defined 
+at $T$-points:
 \begin{equation} \label{Eq_DOM_curl}
 \begin{split}
@@ -110 +157 @@
 \end{equation}
 
-In the special case of pure $z$-coordinates system, \eqref{Eq_DOM_lap} and \eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor becomes a function of the single variable $k$ and thus does not depend on the horizontal location of a grid point. It can be simplified from outside and inside the $\delta _i$ and $\delta_j$ operators. 
-
-The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):
+In the special case of a pure $z$-coordinate system, \eqref{Eq_DOM_lap} and 
+\eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor 
+becomes a function of the single variable $k$ and thus does not depend on the 
+horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to: 
+\begin{equation*}
+\nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} }\left( {\delta 
+_i \left[ {e_{2u} a_1 } \right]+\delta _j \left[ {e_{1v}  a_2 } 
+\right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right]
+\end{equation*}
+
+The vertical average over the whole water column denoted by an overbar becomes 
+for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):
 \begin{equation} \label{DOM_bar}
 \bar q   = \frac{1}{H}\int_{k^b}^{k^o} {q\;e_{3q} \,dk} 
       \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }
 \end{equation}
-where $H_q$  the ocean depth, is the masked sum of the vertical scale factors at q points, $k^b$ and $k^o$ are the bottom and surface $k$-index, and the symbol $k^o$ referring to a summation over all grid points of the same species in the direction indicated by the subscript (here $k$). 
-
-In continuous, the following properties are satisfied:
+where $H_q$  is the ocean depth, which is the masked sum of the vertical scale 
+factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, 
+and the symbol $k^o$ refers to a summation over all grid points of the same type 
+in the direction indicated by the subscript (here $k$). 
+
+In continuous form, the following properties are satisfied:
 \begin{equation} \label{Eq_DOM_curl_grad}
 \nabla \times \nabla q ={\rm {\bf {0}}}
@@ -127 +186 @@
 \end{equation}
 
-It is straight forward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $T$-points and the vector \textbf{A} has its components defined at vector points $(u,v,w)$.
-
-Let $a$ and $b$ be two fields defined on the ocean mesh, extended to zero inside continental area. By integration by part it can be shown that the derivation operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear operators, and further that the averaging operators $\overline{\cdot}^i$, $\overline{\cdot}^j$ and $\overline{\cdot}^k$) are symmetric linear operators, i.e.,
-\begin{equation} \label{DOM_di_adj}
-\sum\limits_i {a_i \;\delta _i \left[ b \right]} \equiv -\sum\limits_i 
-{\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }
-\end{equation}
-\begin{equation} \label{DOM_mi_adj}
-\sum\limits_i {a_i \;\overline b ^i} \equiv \sum\limits_i {\overline a ^{i+1/2}\;b_{i+1/2} } 
-\end{equation}
-
-In other words, the adjoint of the derivation and averaging operators are $\delta_i^*=\delta_{i+1/2}$ and $\overline{\cdot}^{i\,*}= \overline{\cdot}^{i+1/2}$, respectively. These two properties will be used extensively in the \colorbox{yellow} {Appendix C} to 
+It is straightforward to demonstrate that these properties are verified locally in 
+discrete form as soon as the scalar $q$ is taken at $T$-points and the vector 
+\textbf{A} has its components defined at vector points $(u,v,w)$.
+
+Let $a$ and $b$ be two fields defined on the mesh, with value zero inside 
+continental area. Using integration by parts it can be shown that the differencing 
+operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear 
+operators, and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, 
+$\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear 
+operators, $i.e.$
+\begin{align} 
+\label{DOM_di_adj}
+\sum\limits_i { a_i \;\delta _i \left[ b \right]} 
+   &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\
+\label{DOM_mi_adj}
+\sum\limits_i { a_i \;\overline b^{\,i}} 
+   & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 
+\end{align}
+
+In other words, the adjoint of the differencing and averaging operators are 
+$\delta_i^*=\delta_{i+1/2}$ and 
+${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. 
+These two properties will be used extensively in the Appendix~\ref{Apdx_C} to 
 demonstrate integral conservative properties of the discrete formulation chosen.
 
 % -------------------------------------------------------------------------------------------------------------
-%        Numerical Indexation 
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Numerical Indexation}
+%        Numerical Indexing 
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Numerical Indexing}
 \label{DOM_Num_Index}
 
@@ -150 +220 @@
 \begin{figure}[!tb] \label{Fig_index_hor}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./Figures/Fig_index_hor.pdf}
-\caption{Horizontal integer indexation used in the \textsc{Fortran} code. The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices}
+\caption{Horizontal integer indexing used in the \textsc{Fortran} code. The dashed 
+area indicates the cell in which variables contained in arrays have the same 
+$i$- and $j$-indices}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The array representation used in the \textsc{Fortran} code requires an integer indexation while the analytical definition of the mesh (see \S\ref{DOM_cell}) is associated with the use of integer values of $(i,j,k)$ for $T$-points whereas all the other points use both integer and integer and a half values of $(i,j,k)$. Therefore a specific integer indexation must be defined for the latter grid-points 
-(i.e. velocity and vorticity grid-points). Furthermore, it has been chosen to change the direction of the vertical indexation so that the surface level is at $k=1$.
+The array representation used in the \textsc{Fortran} code requires an integer 
+indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is 
+associated with the use of integer values for $T$-points and both integer and 
+integer and a half values for all the other points. Therefore a specific integer 
+indexing must be defined for points other than $T$-points ($i.e.$ velocity and 
+vorticity grid-points). Furthermore, the direction of the vertical indexing has 
+been changed so that the surface level is at $k=1$.
 
 % -----------------------------------
-%        Horizontal Indexation 
+%        Horizontal Indexing 
 % -----------------------------------
-\subsubsection{Horizontal Indexation}
+\subsubsection{Horizontal Indexing}
 \label{DOM_Num_Index_hor}
 
-The indexation in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $T$-point and the eastward $u$-point (northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its nearby northeast $f$-point have the same $i$-and $j$-indices.
+The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $T$-point 
+and the eastward $u$-point (northward $v$-point) have the same index 
+(see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its 
+nearest northeast $f$-point have the same $i$-and $j$-indices.
 
 % -----------------------------------
-%        Vertical Indexation 
+%        Vertical indexing 
 % -----------------------------------
-\subsubsection{Vertical Indexation}
+\subsubsection{Vertical Indexing}
 \label{DOM_Num_Index_vertical}
 
-In the vertical plane, the chosen indexation requires special attention since the $k$-axis is re-oriented downward in the \textsc{Fortran} code compared to the indexation used for the semi-discrete equations and given in \S\ref{DOM_cell}. The sea surface corresponds to the $w$-level $k=1$ like the $T-$level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) is either the ocean bottom or inside the ocean floor while the last $T-$level is always inside the floor (Fig.\ref{Fig_index_vert}). Note that for an increasing $k$ index, a $w$-point and the $T$-point just below have the same $k$ index, in opposition to what is done in the horizontal plane where 
-it is the $T-$point and the nearby 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}). As the scale factors are chosen to be strictly positive, \emph{a minus sign appears in the \textsc{Fortran} code before all the vertical derivatives of the discrete equations given in this documentation}.
+In the vertical, the chosen indexing requires special attention since the 
+$k$-axis is re-orientated downward in the \textsc{Fortran} code compared 
+to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}. 
+The sea surface corresponds to the $w$-level $k=1$ which is the same index 
+as $T$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) 
+either corresponds to the ocean floor or is inside the bathymetry while the last 
+$T$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that 
+for an increasing $k$ index, a $w$-point and the $T$-point just below have the 
+same $k$ index, in opposition to what is done in the horizontal plane where 
+it is the $T$-point and the nearest velocity points in the direction of the horizontal 
+axis that have the same $i$ or $j$ index (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.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!pt] \label{Fig_index_vert}  \begin{center}
 \includegraphics[width=.90\textwidth]{./Figures/Fig_index_vert.pdf}
-\caption{Vertical integer indexation used in the \textsc{Fortran } code. Note that the $k$-axis is oriented downward. The dashed area indicates the cell in which variables contained in arrays have the same $k$-index.}
+\caption{Vertical integer indexing used in the \textsc{Fortran } code. Note that 
+the $k$-axis is orientated downward. The dashed area indicates the cell in 
+which variables contained in arrays have the same $k$-index.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % -----------------------------------
-%        Vertical Indexation 
+%        Domain Size
 % -----------------------------------
-\subsubsection{Domain size}
+\subsubsection{Domain Size}
 \label{DOM_size}
 
-The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpk} in the $i$, $j$ and $k$ directions respectively. They are given as parameters in the \mdl{par\_oce} module (or additional files included in this module such as \textit{par\_ORCA\_R2.h90}, specific to a given configuration). The use of parameters rather than variables (together with dynamic allocation of arrays) was made because it ensured that the compiler would optimize the executable code efficiently, especially on vector machines (optimization may be less efficient when the problem size is unknown at the time of compilation). Nevertheless, it is possible to set up the code with full dynamical allocation by using the AGRIF packaged \colorbox{yellow}{(ref agrif!+ ref part of the doc)}. Note that are other parameters in \mdl{par\_oce} that refer to the domain size. The two parameters $jpidta$ and $jpjdta$, may be larger than $jpiglo$, $jpjglo$ when the user wants to use only a sub-region of a given configuration. This is the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \S\ref{LBC_mpp}).
+The total size of the computational domain is set by the parameters \jp{jpiglo}, 
+\jp{jpjglo} and \jp{jpk} in the $i$, $j$ and $k$ directions respectively. They are 
+given as parameters in the \mdl{par\_oce} module\footnote{When a specific 
+configuration is used (ORCA2 global ocean, etc...) the parameter are actually 
+defined in additional files introduced by \mdl{par\_oce} module via CPP 
+\textit{include} command. For example, ORCA2 parameters are set in 
+\textit{par\_ORCA\_R2.h90} file}. The use of parameters rather than variables 
+(together with dynamic allocation of arrays) was chosen because it ensured that 
+the compiler would optimize the executable code efficiently, especially on vector 
+machines (optimization may be less efficient when the problem size is unknown 
+at the time of compilation). Nevertheless, it is possible to set up the code with full 
+dynamical allocation by using the AGRIF packaged \citep{Debreu_al_CG2008}. 
+%
+\gmcomment{  add the following ref 
+\colorbox{yellow}{(ref part of the doc)} } 
+%
+Note that are other parameters in \mdl{par\_oce} that refer to the domain size. 
+The two parameters $jpidta$ and $jpjdta$ may be larger than $jpiglo$, $jpjglo$ 
+when the user wants to use only a sub-region of a given configuration. This is 
+the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of 
+the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters 
+$jpi$ and $jpj$ refer to the size of each processor subdomain when the code is 
+run in parallel using domain decomposition (\key{mpp\_mpi} defined, see 
+\S\ref{LBC_mpp}).
 
 % ================================================================
 % Domain: Horizontal Grid (mesh) 
 % ================================================================
-\section{Domain: Horizontal Grid (mesh) (\mdl{domhgr} module)}
+\section  [Domain: Horizontal Grid (mesh) (\textit{domhgr})]               
+      {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} }
 \label{DOM_hgr}
 
@@ -201 +319 @@
 \label{DOM_hgr_coord_e}
 
-The ocean mesh (i.e. the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half values of as indicated in Table~\ref{Tab_cell}. The associated scale factors are defined using the analytical first derivative of the transformation \eqref{Eq_scale_factors}. These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which provide the horizontal and vertical meshes, respectively. This section deals with the horizontal mesh parameters.
-
-In a horizontal plane, the location of all the model grid points is defined from the analytical expressions of the latitude $\varphi$ and the longitude $\lambda$ as a function of  $(i,j)$. The horizontal scale factors are calculated using \eqref{Eq_scale_factors}. For example, when the latitude and longitude are function of a single value ($j$ and $i$, respectively) (geographical configuration of the mesh), the horizontal mesh definition reduces to define the wanted $\varphi(j)$, $\varphi'(j)$, $\lambda(i)$, and $\lambda'(i)$ in the \mdl{domhgr} module. The model computes the grid-point positions and scale factors in the horizontal plane as follows:
+The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined 
+by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
+The grid-points are located at integer or integer and a half values of as indicated 
+in Table~\ref{Tab_cell}. The associated scale factors are defined using the 
+analytical first derivative of the transformation \eqref{Eq_scale_factors}. These 
+definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which 
+provide the horizontal and vertical meshes, respectively. This section deals with 
+the horizontal mesh parameters.
+
+In a horizontal plane, the location of all the model grid points is defined from the 
+analytical expressions of the longitude $\lambda$ and  latitude $\varphi$ as a 
+function of  $(i,j)$. The horizontal scale factors are calculated using 
+\eqref{Eq_scale_factors}. For example, when the longitude and latitude are 
+function of a single value ($i$ and $j$, respectively) (geographical configuration 
+of the mesh), the horizontal mesh definition reduces to define the wanted 
+$\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the 
+\mdl{domhgr} module. The model computes the grid-point positions and scale 
+factors in the horizontal plane as follows:
 \begin{flalign*}
-\lambda_T &\equiv \text{glamt} = \lambda(i)  &     \varphi_T &\equiv \text{gphit} = \varphi(j)\\
-\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)&    \varphi_u &\equiv \text{gphiu}= \varphi(j)\\
-\lambda_v &\equiv \text{glamv}= \lambda(i)      &     \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\
-\lambda_f &\equiv \text{glamf }= \lambda(i+1/2)&      \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) 
+\lambda_T &\equiv \text{glamt}= \lambda(i)     & \varphi_T &\equiv \text{gphit} = \varphi(j)\\
+\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\
+\lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\
+\lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) 
 \end{flalign*}
 \begin{flalign*}
 e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |&
-e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)|\\
+e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\
 e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |&
 e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\
@@ -220 +353 @@
 e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|
 \end{flalign*}
-where the last letter of each computational name indicates the grid point considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). Note that the horizontal position and scale factors of $w$-points are exactly equal to those of $T-$points, thus no specific arrays are defined at those grid-points. 
-
-Note that the definition of the scale factors --- as the analytical first derivative of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$ --- is specific to the OPA model \citep{Marti1992}. As an example, $e_{1T}$ is defined locally at a $T$-point, whereas many other models on a C grid choose to define such a scale factor as the distance between the $U$-points on each side of the $T$-point. Relying on an analytical transformation has two advantages: firstly, there is no ambiguity in the scale factors appearing in the discrete equations, since they are first introduced in the continuous equations; secondly, analytical transformations encourage good practice by the definition of smooth grids \citep{Treguier1996}. An example of the effect of such a choice is shown in Fig.~\ref{Fig_zgr_e3}.
+where the last letter of each computational name indicates the grid point 
+considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with 
+all universal constants). Note that the horizontal position of and scale factors 
+at $w$-points are exactly equal to those of $T$-points, thus no specific arrays 
+are defined at $w$-points. 
+
+Note that the definition of the scale factors ($i.e.$ as the analytical first derivative 
+of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is 
+specific to the \NEMO model \citep{Marti1992}. As an example, $e_{1T}$ is defined 
+locally at a $T$-point, whereas many other models on a C grid choose to define 
+such a scale factor as the distance between the $U$-points on each side of the 
+$T$-point. Relying on an analytical transformation has two advantages: firstly, there 
+is no ambiguity in the scale factors appearing in the discrete equations, since they 
+are first introduced in the continuous equations; secondly, analytical transformations 
+encourage good practice by the definition of smoothly varying grids (rather than 
+allowing the user to set arbitrary jumps in thickness between adjacent layers) 
+\citep{Treguier1996}. An example of the effect of such a choice is shown in 
+Fig.~\ref{Fig_zgr_e3}.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_zgr_e3}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./Figures/Fig_zgr_e3.pdf}
-\caption{(a) Traditional definition of grid-point position and grid-size in the vertical versus (b) analytically derived grid-point position and scale factors. For both grid here,a same $w$-point depth has been chosen but in (a) the $T$-points are set at the middle of $w$-points while in (b) they are defined from an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. Note the resulting difference between the value of the grid-size $\Delta_k$ and those of the scale factor $e_k$. }
+\caption{Comparison of (a) traditional definitions of grid-point position and grid-size 
+in the vertical, and (b) analytically derived grid-point position and scale factors. For 
+both grids here,  the same $w$-point depth has been chosen but in (a) the 
+$T$-points are set half way between $w$-points while in (b) they are defined from 
+an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. 
+Note the resulting difference between the value of the grid-size $\Delta_k$ and 
+those of the scale factor $e_k$. }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % -------------------------------------------------------------------------------------------------------------
-%        Choice of horizontal grids 
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Choice of horizontal grids}
+%        Choice of horizontal grid
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Choice of horizontal grid}
 \label{DOM_hgr_msh_choice}
 
-The user has three options to define a horizontal grid, involving the parameter $jphgr\_mesh$ of the \mdl{par\_oce} module. 
-\begin{enumerate}
-\item For the most general curvilinear orthogonal grids, the coordinates and their first derivatives with respect to $i$ and $j$ are provided in a file, read in \rou{hgr\_read} subroutine of the domhgr module: \jp{jphgr\_mesh}=0.
-\item A few simple analytical grids are provided as examples, that can be selected by setting \jp{jphgr\_mesh}=1 to 5 (see below) 
-\item For other analytical grids, the \mdl{domhgr} module must be modified by the user. 
-\end{enumerate}
-
-There are two simple cases of geographical grids on the sphere. With \jp{jphgr\_mesh}=1, the grid is regular in space, with grid sizes specified by parameters \pp{ppe1\_deg} and \pp{ppe2\_deg}, respectively. A geographical grid 
-can be very anisotropic at high latitudes, because of the convergence of meridians (the zonal scale factors $e_1$ become much smaller than the meridional scale factors $e_2$). The Mercator grid (\jp{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale factors in the same way as the zonal ones. In that case, meridional scale factors and latitudes are calculated analytically using the formulae appropriate for a Mercator projection, based on \pp{ppe1\_deg} which is a reference grid spacing at the equator (this applies even when the geographical equator is situated outside the model domain). In those two cases (\jp{jphgr\_mesh}=1 or 4), the grid position is defined by the longitude and latitude of the south-westhernmost point (\pp{ppglamt0} and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide an approximate starting latitude: the real latitude will be recalculated analytically, so as to ensure that the equator corresponds to a $T$- and$ U$-point.  
-
-Rectangular grids ignoring the spherical geometry are defined with \jp{jphgr\_mesh} = 2, 3, 5. The domain is either a $f$-plane (\jp{jphgr\_mesh} = 2, Coriolis factor is constant) or a beta-plane (\jp{jphgr\_mesh} = 3, the Coriolis factor is linear in the $j$-direction). The grid size is uniform in each direction, and given in meters by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively. The zonal grid coordinate (glam. arrays) is in kilometers, starting at zero with the first T point. The meridional coordinate (gphi. arrays) is in kilometers, and the second $T$-point corresponds to coordinate gphit=0. The input parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference latitude for computation of the Coriolis parameter. In the case of the beta plane, \pp{ppgphi0} corresponds to the center of the domain. Finally, the special case \jp{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the GYRE configuration representing a classical mid-latitude double gyre system. The rotation allows to maximize the jet length relative to the gyre areas (and the number of grid points). 
-
-The choice of the grid must be consistent with the boundary conditions specified by the parameter \jp{jperio} (see {\S\ref{LBC}).
+The user has three options available in defining a horizontal grid, which involve 
+the parameter $jphgr\_mesh$ of the \mdl{par\_oce} module. 
+\begin{description}
+\item[\jp{jphgr\_mesh}=0]  The most general curvilinear orthogonal grids.
+The coordinates and their first derivatives with respect to $i$ and $j$ are 
+provided in a file, read in \rou{hgr\_read} subroutine of the domhgr module.
+\item[\jp{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below). 
+For other analytical grids, the \mdl{domhgr} module must be modified by the user. 
+\end{description}
+
+There are two simple cases of geographical grids on the sphere. With 
+\jp{jphgr\_mesh}=1, the grid (expressed in degrees) is regular in space, 
+with grid sizes specified by parameters \pp{ppe1\_deg} and \pp{ppe2\_deg}, 
+respectively. Such a geographical grid can be very anisotropic at high latitudes 
+because of the convergence of meridians (the zonal scale factors $e_1$ 
+become much smaller than the meridional scale factors $e_2$). The Mercator 
+grid (\jp{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale 
+factors in the same way as the zonal ones. In this case, meridional scale factors 
+and latitudes are calculated analytically using the formulae appropriate for 
+a Mercator projection, based on \pp{ppe1\_deg} which is a reference grid spacing 
+at the equator (this applies even when the geographical equator is situated outside 
+the model domain). 
+%%%
+\gmcomment{ give here the analytical expression of the Mercator mesh}
+%%%
+In these two cases (\jp{jphgr\_mesh}=1 or 4), the grid position is defined by the 
+longitude and latitude of the south-westernmost point (\pp{ppglamt0} 
+and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide 
+an approximate starting latitude: the real latitude will be recalculated analytically, 
+in order to ensure that the equator corresponds to line passing through $T$- 
+and $u$-points.  
+
+Rectangular grids ignoring the spherical geometry are defined with 
+\jp{jphgr\_mesh} = 2, 3, 5. The domain is either an $f$-plane (\jp{jphgr\_mesh} = 2, 
+Coriolis factor is constant) or a beta-plane (\jp{jphgr\_mesh} = 3, the Coriolis factor 
+is linear in the $j$-direction). The grid size is uniform in meter in each direction, 
+and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively. 
+The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero 
+with the first $T$-point. The meridional coordinate (gphi. arrays) is in kilometers, 
+and the second $T$-point corresponds to coordinate $gphit=0$. The input 
+parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference 
+latitude for computation of the Coriolis parameter. In the case of the beta plane, 
+\pp{ppgphi0} corresponds to the center of the domain. Finally, the special case 
+\jp{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the 
+GYRE configuration, representing a classical mid-latitude double gyre system. 
+The rotation allows us to maximize the jet length relative to the gyre areas 
+(and the number of grid points). 
+
+The choice of the grid must be consistent with the boundary conditions specified 
+by the parameter \jp{jperio} (see {\S\ref{LBC}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -256 +447 @@
 \label{DOM_hgr_files}
 
-All the arrays related to a particular ocean model configuration (grid-point position, scale factors, masks) can be saved in files if $\np{nmsh} \not= 0$ (namelist parameter). This can be particularly useful for plots and off-line diagnostics. In some cases, the user may choose to make a local modification of a scale factor in the code. This is the case in global configurations when restricting the width of a specific strait (usually a one-grid-point strait that happens to be too wide due to the insufficient model resolution). On example is Lombok Strait in the ORCA2 configuration. When such modifications are done, the output grid written when $\np{nmsh} \not=0$ is not exactly equal to the input grid.
+All the arrays relating to a particular ocean model configuration (grid-point 
+position, scale factors, masks) can be saved in files if $\np{nmsh} \not= 0$ 
+(namelist parameter). This can be particularly useful for plots and off-line 
+diagnostics. In some cases, the user may choose to make a local modification 
+of a scale factor in the code. This is the case in global configurations when 
+restricting the width of a specific strait (usually a one-grid-point strait that 
+happens to be too wide due to insufficient model resolution). An example 
+is Gibraltar Strait in the ORCA2 configuration. When such modifications are done, 
+the output grid written when $\np{nmsh} \not=0$ is no more equal to the input grid.
 
 % ================================================================
 % Domain: Vertical Grid (domzgr)
 % ================================================================
-\section{Domain: Vertical Grid (\mdl{domzgr} module)}
+\section  [Domain: Vertical Grid (\textit{domzgr})]
+      {Domain: Vertical Grid \small{(\mdl{domzgr} module)} }
 \label{DOM_zgr}
 %-----------------------------------------nam_zgr & namdom-------------------------------------------
@@ -268 +468 @@
 %-------------------------------------------------------------------------------------------------------------
 
-In the vertical, the model mesh is determined by four things: (1) the bathymetry given in meters ; (2) the number of levels of the model (\jp{jpk}) ; (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors 
-(derivatives of the transformation) ; and (4) the masking system, i.e. the number of wet model levels at each $(i,j)$.
+In the vertical, the model mesh is determined by four things: 
+(1) the bathymetry given in meters ; 
+(2) the number of levels of the model (\jp{jpk}) ; 
+(3) the analytical transformation $z(i,j,k)$ and the vertical scale factors 
+(derivatives of the transformation) ; 
+and (4) the masking system, $i.e.$ the number of wet model levels at each 
+$(i,j)$ column of points.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!tb] \label{Fig_z_zps_s_sps}  \begin{center}
 \includegraphics[width=1.0\textwidth]{./Figures/Fig_z_zps_s_sps.pdf}
-\caption{The ocean bottom as seen by the model: (a) $z$-coordinate with full step, (b) $z$-coordinate with partial step, (c) $s$-coordinate: terrain following representation, (d) hybrid $s-z$ coordinate, (e) hybrid $s-z$ coordinate with partial step, and (f) same as (e) but with variable volume level associated with the non-linear free surface. Note that the variable volume level (\key{vvl}) could be used with any of the 5 coordinates (a) to (e).}
+\caption{The ocean bottom as seen by the model: 
+(a) $z$-coordinate with full step, 
+(b) $z$-coordinate with partial step, 
+(c) $s$-coordinate: terrain following representation, 
+(d) hybrid $s-z$ coordinate, 
+(e) hybrid $s-z$ coordinate with partial step, and 
+(f) same as (e) but with variable volume associated with the non-linear free surface. 
+Note that the variable volume option (\key{vvl}) can be used with any of the 
+5 coordinates (a) to (e).}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The choice of a vertical coordinate among all those offered in NEMO, even if it is made through a namelist parameter, must be done once of all at the beginning of an experiment. It is not intended as an option which can be 
-enabled or disabled in the middle of an experiment. Three main choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step bathymetry (\np{ln\_zco}=true), $z$-coordinate with partial step bathymetry (\np{ln\_zps}=true), or generalized, $s$-coordinate (\np{ln\_sco}=true). Hybridation of the three main 
-coordinates are available: hybrid $s-z$ or $s-zps$ coordinate (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable volume option \key{vvl}), the coordinate follow the time-variation of the free surface so that the transformation is time dependent: $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step bathymetry or $s$-coordinates (hybride and partial step coordinates not yet implemented in NEMO v2.3). 
-
-Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for reading it from a file. The only input file is the bathymetry (in meters). After reading the bathymetry, the algorithm for vertical grid definition differs between the different options:
+The choice of a vertical coordinate, even if it is made through a namelist parameter, 
+must be done once of all at the beginning of an experiment. It is not intended as an 
+option which can be enabled or disabled in the middle of an experiment. Three main 
+choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step 
+bathymetry (\np{ln\_zco}=true), $z$-coordinate with partial step bathymetry 
+(\np{ln\_zps}=true), or generalized, $s$-coordinate (\np{ln\_sco}=true). 
+Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate 
+(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable 
+volume option \key{vvl}) ($i.e.$ non-linear free surface), the coordinate follow the 
+time-variation of the free surface so that the transformation is time dependent: 
+$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step 
+bathymetry or $s$-coordinate (hybride and partial step coordinates have not 
+yet been tested in NEMO v2.3). 
+
+Contrary to the horizontal grid, the vertical grid is computed in the code and no 
+provision is made for reading it from a file. The only input file is the bathymetry 
+(in meters)\footnote{N.B. in full step $z$-coordinate, a \textit{bathy\_level} file can 
+replace the \textit{bathy\_meter} file, so that the computation of the number of 
+wet ocean point in each water column is by-passed}. After reading the bathymetry, 
+the algorithm for vertical grid definition differs between the different options:
 \begin{description}
 \item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$.
 \item[\textit{zps}] set a reference coordinate transformation $z_0 (k)$, and 
-calculate the height at the deepest levels using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
-\item[\textit{sco}] Smooth the bathymetry to fullfill the hydrostatic consistency criteria and set the three-dimensional transformation.
-\item[\textit{s-z} and \textit{s-zps}] Smooth the bathymetry to fullfill the hydrostatic consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and possibly introduce masking of extra land points to better fit the original bathymetry file
+calculate the thickness of the deepest level at each $(i,j)$ point using the 
+bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
+\item[\textit{sco}] smooth the bathymetry to fulfil the hydrostatic consistency 
+criteria and set the three-dimensional transformation.
+\item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfil the hydrostatic 
+consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and 
+possibly introduce masking of extra land points to better fit the original bathymetry file
 \end{description}
-
-Generally, the arrays describing the grid point depths and vertical scale factors are three dimensional arrays $(i,j,k)$. In the special case of $z$-coordinates with full step bottom topography, it is possible to define 
-those arrays as one-dimensional, in order to save memory. This is performed by defining the \key{zco} C-Pre-Processor (CPP) key. To improve the code readability while providing this flexibility, the vertical coordinates and scale factors are defined as functions of $(i,j,k)$ with "fs" as prefix (examples: \textit{fsdeptht, fse3t,} etc) that can be equal to three-dimensional arrays, or a one dimensional array when \key{zco} is defined. These functions are defined in the file 
-\textit{domzgr\_substitute.h90} of the DOM directory. They are used through the code, and replaced by the corresponding arrays at the time of pre-processing (CPP capability).
+%%%
+\gmcomment{   add the description of the smoothing:  envelop topography...}
+%%%
+
+Generally, the arrays describing the grid point depths and vertical scale factors 
+are three dimensional arrays $(i,j,k)$. In the special case of $z$-coordinates with 
+full step bottom topography, it is possible to define those arrays as one-dimensional, 
+in order to save memory. This is performed by defining the \key{zco} 
+C-Pre-Processor (CPP) key. To improve the code readability while providing this 
+flexibility, the vertical coordinate and scale factors are defined as functions of 
+$(i,j,k)$ with "fs" as prefix (examples: \textit{fsdeptht, fse3t,} etc) that can be either 
+three-dimensional arrays, or a one dimensional array when \key{zco} is defined. 
+These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory. 
+They are used throughout the code, and replaced by the corresponding arrays at 
+the time of pre-processing (CPP capability).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -304 +547 @@
 namelist variable \np{ntopo}: 
 \begin{description}
-\item[\np{ntopo} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ is given by the coordinate transformation. The domain can either a closed basin or a periodic channel according to the parameter \jp{jperio}. 
-\item[\np{ntopo} = -1] a domain with a bump of topography at the central latitude and 1/3 of the domain width. This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 
-\item[\np{ntopo} = 1] read a bathymetry. The bathymetry file (Netcdf format) provides the ocean depth (positive, in meters) at each grid point of the model grid. The bathymetry is usually built by interpolating a standard bathymetry product (e.g., ETOPO2) onto the horizontal ocean mesh. The bathymetry file defines the coastline: where the bathymetry is zero, no model levels are defined (all levels are masked).
+\item[\np{ntopo} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ 
+is given by the coordinate transformation. The domain can either be a closed 
+basin or a periodic channel depending on the parameter \jp{jperio}. 
+\item[\np{ntopo} = -1] a domain with a bump of topography one third of the 
+domain width at the central latitude. This is meant for the "EEL-R5" configuration, 
+a periodic or open boundary channel with a seamount. 
+\item[\np{ntopo} = 1] read a bathymetry. The bathymetry file (Netcdf format) 
+provides the ocean depth (positive, in meters) at each grid point of the model grid. 
+The bathymetry is usually built by interpolating a standard bathymetry product 
+($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also 
+defines the coastline: where the bathymetry is zero, no model levels are defined 
+(all levels are masked).
 \end{description}
 
-When using the rigid lid approximation (\key{dynspg\_rl} defined) isolated land masses (islands) must be identified by negative integers in the input bathymetry file (see \S\ref{MISC_solisl}).
-
-When the ocean is coupled to an atmospheric model it is better to represent 
+When using the rigid lid approximation (\key{dynspg\_rl} is defined) isolated land 
+masses (islands) must be identified by negative integers in the input bathymetry file 
+(see \S\ref{MISC_solisl}).
+
+When a global ocean is coupled to an atmospheric model it is better to represent 
 all large water bodies (e.g, great lakes, Caspian sea...) even if the model 
-resolution does not allow to represent their communication with the rest of 
-the ocean. This is unnecessary when the ocean is forced by fixed atmospheric 
-conditions. A possibility is offered to the user to set to zero the 
-bathymetry in rectangular regions covering those closed seas (see \S\ref{MISC_closea}), but the code has to be adapted to the user's configuration. 
+resolution does not allow their communication with the rest of the ocean. 
+This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
+so these seas can be removed from the ocean domain. The user has the option 
+to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the 
+code has to be adapted to the user's configuration. 
 
 % -------------------------------------------------------------------------------------------------------------
 %        z-coordinate  and reference coordinate transformation
 % -------------------------------------------------------------------------------------------------------------
-\subsection [$z$-coordinate (\np{ln\_zco}=T or \key{zco})]
-         {$z$-coordinate (\np{ln\_zco}=T or \key{zco}) and reference coordinate}
+\subsection[$z$-coordinate (\np{ln\_zco} or \key{zco})]
+        {$z$-coordinate (\np{ln\_zco}=.true. or \key{zco}) and reference coordinate}
 \label{DOM_zco}
 
@@ -328 +583 @@
 \begin{figure}[!tb] \label{Fig_zgr}  \begin{center}
 \includegraphics[width=0.90\textwidth]{./Figures/Fig_zgr.pdf}
-\caption{Default vertical mesh for ORCA2-L30. Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from (III.2.1) in z-coordinates.}
+\caption{Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level 
+functions for (a) T-point depth and (b) the associated scale factor as computed 
+from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The reference coordinate transformation $z_0 (k)$ defines the arrays \textit{gdept0} and \textit{gdepw0} for $T$- and $w$-points, respectively. As indicated on Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw(1)$ being the ocean surface. There are at most \jp{jpk}-1 $T$-points in the ocean, the additional $T$-point at $jk=jpk$ is below the sea floor and is not used. The vertical location of $w$- and $T$-levels is defined from the analytic expression of the depth $z_0 (k)$ whose analytic derivative with respect to $k$ provides the vertical scale factors. The user must provide the analytical expression of both $z_0 $and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} through statement functions, using parameters provided in the \textit{par\_oce.h90} file. 
+The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ 
+and $gdepw_0$ for $T$- and $w$-points, respectively. As indicated on 
+Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the 
+ocean surface. There are at most \jp{jpk}-1 $T$-points inside the ocean, the 
+additional $T$-point at $jk=jpk$ is below the sea floor and is not used. 
+The vertical location of $w$- and $T$-levels is defined from the analytic expression 
+of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the 
+vertical scale factors. The user must provide the analytical expression of both 
+$z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} 
+through statement functions, using parameters provided in the \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. 
 
-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 $z$-coordinates and partial steps: 
+For climate-related studies it is often desirable to concentrate the vertical resolution 
+near the ocean surface. The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 
 \begin{equation} \label{DOM_zgr_ana}
 \begin{split}
@@ -343 +610 @@
 \end{split}
 \end{equation}
-where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with a smooth hyperbolic tangent transition in between (Fig.~\ref{Fig_zgr}).
-
-The first grid defined for ORCA2 had $10~m$ ($500~m)$ resolution in the surface (bottom) layers and a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. This leads to the following conditions:
+where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an 
+expression allows us to define a nearly uniform vertical location of levels at the 
+ocean top and bottom with a smooth hyperbolic tangent transition in between 
+(Fig.~\ref{Fig_zgr}).
+
+The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the 
+surface (bottom) layers and a depth which varies from 0 at the sea surface to a 
+minimum of $-5000~m$. This leads to the following conditions:
 \begin{equation} \label{DOM_zgr_coef}
 \begin{split}
@@ -355 +627 @@
 \end{equation}
 
-With the choice of the stretching $h_{cr} =3$ and the number of levels \jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in \eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. For the first standard 
-ORCA2 vertial grid this led to the following values: $h_{sur} =4762.96$, $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters \pp{ppsur}, \pp{ppa0}, \pp{ppa1}, \pp{ppkth} in the parameter file \mdl{par\_oce}. 
+With the choice of the stretching $h_{cr} =3$ and the number of levels 
+\jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in 
+\eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is 
+satisfied, through an optimisation procedure using a bisection method. For the first 
+standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$, 
+$h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and 
+scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and 
+given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters 
+\pp{ppsur}, \pp{ppa0}, \pp{ppa1}, \pp{ppkth} in the parameter file \mdl{par\_oce}. 
 
 Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is 
 possible to recalculate them. In that case the user sets 
-\pp{ppsur}=\pp{ppa0}=\pp{ppa1}=\pp{pp\_to\_be\_computed}, in \mdl{par\_oce}, and specifies instead the four following parameters:
+\pp{ppsur}=\pp{ppa0}=\pp{ppa1}=\pp{pp\_to\_be\_computed}, in \mdl{par\_oce}, 
+and specifies instead the four following parameters:
 \begin{itemize}
-\item    \pp{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger \pp{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual.
-\item    \pp{ppkth}=$h_{th} $: is approximately the model level at which maximum stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
+\item    \pp{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger 
+\pp{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual.
+\item    \pp{ppkth}=$h_{th} $: is approximately the model level at which maximum 
+stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
 \item    \pp{ppdzmin}: minimum thickness for the top layer (in meters)
 \item    \pp{pphmax}: total depth of the ocean (meters).
 \end{itemize}
-As an example, for the $45$ layers used in DRAKKAR configuration those parameters are: \jp{jpk}=46, \pp{ppacr}=9, \pp{ppkth}=23.563, \pp{ppdzmin}=6m, \pp{pphmax}=5750m.
+As an example, for the $45$ layers used in the DRAKKAR configuration those 
+parameters are: \jp{jpk}=46, \pp{ppacr}=9, \pp{ppkth}=23.563, \pp{ppdzmin}=6m, 
+\pp{pphmax}=5750m.
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -396 +680 @@
 21 &  \textbf{732.20}   &   611.89 &   \textbf{261.03} &    220.35 \\ \hline
 22 &  \textbf{1033.22}&  872.87 &   \textbf{339.39} &    301.42 \\ \hline
-23 &  \textbf{1405.70}& 1211.59 &   \textbf{402.26} &    373.31 \\ \hline
-24 &  \textbf{1830.89}& 1612.98 &   \textbf{444.87} &    426.00 \\ \hline
-25 &  \textbf{2289.77}& 2057.13 &   \textbf{470.55} &    459.47 \\ \hline
-26 &  \textbf{2768.24}& 2527.22 &   \textbf{484.95} &    478.83 \\ \hline
-27 &  \textbf{3257.48}& 3011.90 &   \textbf{492.70} &    489.44 \\ \hline
-28 &  \textbf{3752.44}& 3504.46 &   \textbf{496.78} &    495.07 \\ \hline
-29 &  \textbf{4250.40}& 4001.16 &   \textbf{498.90} &    498.02 \\ \hline
-30 &  \textbf{4749.91}& 4500.02 &   \textbf{500.00} & 499.54 \\ \hline
+23 &  \textbf{1405.70}& 1211.59 & \textbf{402.26} &   373.31 \\ \hline
+24 &  \textbf{1830.89}& 1612.98 & \textbf{444.87} &   426.00 \\ \hline
+25 &  \textbf{2289.77}& 2057.13 & \textbf{470.55} &   459.47 \\ \hline
+26 &  \textbf{2768.24}& 2527.22 & \textbf{484.95} &   478.83 \\ \hline
+27 &  \textbf{3257.48}& 3011.90 & \textbf{492.70} &   489.44 \\ \hline
+28 &  \textbf{3752.44}& 3504.46 & \textbf{496.78} &   495.07 \\ \hline
+29 &  \textbf{4250.40}& 4001.16 & \textbf{498.90} &   498.02 \\ \hline
+30 &  \textbf{4749.91}& 4500.02 & \textbf{500.00} &   499.54 \\ \hline
 31 &  \textbf{5250.23}& 5000.00 &   \textbf{500.56} & 500.33 \\ \hline
 \end{tabular} \end{center} 
-\caption{Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}}
+\caption{Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration 
+as computed from \eqref{DOM_zgr_ana} using the coefficients given in 
+\eqref{DOM_zgr_coef}}
 \end{table}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -413 +699 @@
 %        z-coordinate with partial step
 % -------------------------------------------------------------------------------------------------------------
-\subsection{z-coordinate with partial step (\np{ln\_zps}=T)}
+\subsection   [$z$-coordinate with partial step (\np{ln\_zps})]
+         {$z$-coordinate with partial step (\np{ln\_zps}=.true.)}
 \label{DOM_zps}
-%--------------------------------------------namdom-----------------------------------------------------
+%--------------------------------------------namdom-------------------------------------------------------
 \namdisplay{namdom} 
 %--------------------------------------------------------------------------------------------------------------
 
-In that case, the depths of the model levels are still defined by the 
+In $z$-coordinate partial step, the depths of the model levels are defined by the 
 reference analytical function $z_0 (k)$ as described in the previous 
-section, excepted in the bottom layer. The thickness of the bottom layer is 
+section, \emph{except} in the bottom layer. The thickness of the bottom layer is 
 allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a 
 better representation of the bathymetry, especially in the case of small 
@@ -431 +718 @@
 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 layers grid, the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). Two 
+\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 
@@ -443 +730 @@
 %        s-coordinate
 % -------------------------------------------------------------------------------------------------------------
-\subsection{$s$-coordinate (\np{ln\_sco}=T)}
+\subsection   [$s$-coordinate (\np{ln\_sco})]
+         {$s$-coordinate (\np{ln\_sco}=true)}
 \label{DOM_sco}
-%--------------------------------------------nam_zgr_sco---------------------------------------------------
+%------------------------------------------nam_zgr_sco---------------------------------------------------
 \namdisplay{nam_zgr_sco} 
 %--------------------------------------------------------------------------------------------------------------
-In s-coordinate (\key{sco} defined), the depths of the model 
-levels are defined from the product of a depth field and a stretching 
-function and its derivative, respectively:
+In $s$-coordinate (\key{sco} is defined), the depth and thickness of the model 
+levels are defined from the product of a depth field and either a stretching 
+function or its derivative, respectively:
 \begin{equation} \label{DOM_sco_ana}
 \begin{split}
- z(k)    &= h(i,j) \; z_0(k)  \\
+ z(k)       &= h(i,j) \; z_0(k)  \\
  e_3(k)  &= h(i,j) \; z_0'(k)
 \end{split}
 \end{equation}
-where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at $T-$point location 
-in the horizontal and $z_0 (k)$ is a function which varies from $0$ at the sea 
+where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $T$-point 
+location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea 
 surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean 
-depth as a mixed step-like and bottom following representation of the 
+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.
 
@@ -472 +760 @@
 \end{split}
 \end{equation}
-where $h_c $is the thermocline depth and $\theta $ and $b$ are the surface and 
+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$. Examples of the stretching function applied to a seamount are given in 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]{./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 resolution}
+\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}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -489 +777 @@
 \label{DOM_zgr_vvl}
 
-This option is described in the report by Levier \textit{et al.} (2007), available on the NEMO web site. 
+This option is described in the Report by Levier \textit{et al.} (2007), available on 
+the \NEMO web site. 
 
 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
@@ -511 +800 @@
 
 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} 
-routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points that do not 
-communicate with another ocean point at the same level are eliminated.
-
-In case of rigid-lid approximation and islands in the computational domain (\np{ln\_dynspg\_rl}=true and \key{island} defined), the \textit{mbathy} array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the 
-following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T-$points are land points of the $n^{th}$ island ; $mbathy(i,j) =0$, $T-$points are land points of the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points are ocean points, the others land points. This is used to compute the island barotropic stream function used in rigid lid computation (see \S\ref{MISC_solisl}).
+routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points 
+that do not communicate with another ocean point at the same level are eliminated.
+
+In the case of the rigid-lid approximation when islands occur in the computational 
+domain (\np{ln\_dynspg\_rl}=.true. and \key{island} is defined), the \textit{mbathy} 
+array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the 
+following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T$-points are 
+land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $T$-points are land points 
+on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points 
+are ocean points, the others are points below the ocean floor. 
+
+This is used to compute the island barotropic stream function used in the rigid lid 
+computation (see \S\ref{MISC_solisl}).
 
 From the \textit{mbathy} array, the mask fields are defined as follows:
 \begin{align*}
-tmask(i,j,k) &= \begin{cases}   1&   \text{ if $k\leq mbathy(i,j)$  }    \\
-                                              0&   \text{ if $k\leq mbathy(i,j)$  }    \end{cases}     \\
-umask(i,j,k) &= \; tmask(i,j,k) \;.\; tmask(i+1,j,k)  \\
-umask(i,j,k) &= \; tmask(i,j,k) \;.\; tmask(i,j+1,k)  \\
-umask(i,j,k) &= \; tmask(i,j,k) \;.\; tmask(i+1,j,k)  \\
-                  & \quad . \; tmask(i,j,k) \;.\; tmask(i+1,j,k)
+tmask(i,j,k) &= \begin{cases}   \; 1&   \text{ if $k\leq mbathy(i,j)$  }    \\
+                                                \; 0&   \text{ if $k\leq mbathy(i,j)$  }    \end{cases}     \\
+umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
+umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\
+umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
+                   & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k)
 \end{align*}
 
-Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with the numerical 
-indexation used (\S~\ref{DOM_Num_Index}). Moreover, the specification of closed lateral 
-boundaries requires that at least the first and last rows and columns of 
-\textit{mbathy} array are set to zero. In the particular case of an east-west cyclic 
-boundary condition, \textit{mbathy} has its last column equal to the second one and its 
-first column equal to the last but one (and so the mask arrays) (see \S~\ref{LBC_jperio}).
-
-\colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.
+\gmcomment{ STEVEN: are the dots multiplications?}      
+
+Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with 
+the numerical indexing used (\S~\ref{DOM_Num_Index}). Moreover, the 
+specification of closed lateral boundaries requires that at least the first and last 
+rows and columns of the \textit{mbathy} array are set to zero. In the particular 
+case of an east-west cyclical boundary condition, \textit{mbathy} has its last 
+column equal to the second one and its first column equal to the last but one 
+(and so too the mask arrays) (see \S~\ref{LBC_jperio}).
+
+%%%
+\gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  }
+%%%
 
 % ================================================================
-% Time discretisation
+% Time Discretisation
 % ================================================================
 \section{Time Discretisation}
 \label{DOM_nxt}
 
-The time stepping used in OPA is a three level scheme that can be presented as follows:
+The time stepping used in \NEMO is a three level scheme that can be 
+represented as follows:
 \begin{equation} \label{Eq_DOM_nxt}
    x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t-\Delta t,t,t+\Delta t}
 \end{equation} 
-where $x$ stand for $u$, $v$, $T$ or $S$, RHS is the Right-Hand-Side of the corresponding time 
-evolution equation, $\Delta t$ is the time step and the overscripts indicate 
-the time at which a quantity is evaluated. Each term of the RHS is evaluated at 
-specific time step(s) depending on the physics to which it is associated. 
+where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the 
+corresponding time evolution equation; $\Delta t$ is the time step; and the 
+superscripts indicate the time at which a quantity is evaluated. Each term of the 
+RHS is evaluated at a specific time step(s) depending on the physics with which 
+it is associated. 
+
 The choice of the time step used for this evaluation is discussed below as 
-well as the implication in term of starting or restarting a model 
+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 successively.
-
-The three level scheme requires three arrays for the prognostic variables. For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array, although referred to as $x_a$ (after) in the code, is usually not the variable $x_a$ at the next time step; rather, it is used to store the time derivative (RHS in \eqref{Eq_DOM_nxt}) prior to time-stepping the equation. Generally, the time stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt} modules, excepted for implicit vertical diffusion or sea surface height when time-splitting options are used.
+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?}
+%%%
+
+The three level scheme requires three arrays for each prognostic variables. 
+For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array, 
+although referred to as $x_a$ (after) in the code, is usually not the variable at 
+the next time step; but rather it is used to store the time derivative (RHS in 
+\eqref{Eq_DOM_nxt}) prior to time-stepping the equation. Generally, the time 
+stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt} 
+modules, except for implicit vertical diffusion or sea surface height when 
+time-splitting options are used.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -564 +879 @@
 
 The time stepping used for non-diffusive processes is the well-known 
-leapfrog scheme. It is a time centred scheme, i.e. the RHS are evaluated at 
+leapfrog scheme. It is a time centred scheme, i.e. the RHS is evaluated at 
 time step $t$, the now time step. It is only used for non-diffusive terms, 
-that is momentum and tracer advection, pressure gradient, and coriolis 
+that is momentum and tracer advection, pressure gradient, and Coriolis 
 terms. This scheme is widely used for advective processes in low-viscosity 
 fluids. It is an efficient method that achieves second-order accuracy with 
@@ -578 +893 @@
 "time-splitting" \citep{Haltiner1980} that develops when the method 
 is used to model non linear fluid dynamics: the even and odd time steps tend 
-to diverge between a physical and a computational mode. Time splitting can 
+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 
+\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 OPA we follow the first 
+integration step with a two-level scheme. In \NEMO we follow the first 
 strategy:
 \begin{equation} \label{Eq_DOM_nxt_asselin}
 x_F^t  = x^t + \gamma \, \left[ x_f^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right]
 \end{equation} 
-where the subscript $f$ denotes filtered values and $\gamma$ is the asselin coefficient. $\gamma$ is initialized as \np{atfp} (namelist parameter). Its default value is \np{atfp}=0.1.  This default value causes a significant dissipation of high frequency motions. Recommanded values in idealized studies of shallow water turbulence are two order of magnitude lower (\citep{Farge1987}). Both strategies do, nevertheless, degrade the accuracy of the calculation from second to first order. The leapfrog scheme associated to a Robert-Asselin 
+where the subscript $f$ denotes filtered values and $\gamma$ is the Asselin 
+coefficient. $\gamma$ is initialized as \np{atfp} (namelist parameter). 
+Its default value is \np{atfp}=0.1.  This default value causes a significant dissipation 
+of high frequency motions. Recommended values in idealized studies of shallow 
+water turbulence are two orders of magnitude smaller (\citep{Farge1987}). 
+Both strategies do, nevertheless, degrade the accuracy of the calculation from 
+second to first order. The leapfrog scheme combined with a Robert-Asselin 
 time filter has been preferred to other time differencing schemes such as 
-predictor corrector or trapezoidal schemes because the user can better 
-control the magnitude and the spatial structure of the time diffusion of the 
-scheme. In association with the centred space discretisation of the 
+predictor corrector or trapezoidal schemes, because the user has an explicit 
+and simple control of the magnitude of the time diffusion of the scheme. 
+In association with the 2nd order centred space discretisation of the 
 advective terms in the momentum and tracer equations, it avoids implicit 
-numerical diffusion in both the time and space discretisation of the 
-advective term: they are both set explicitly by the user through the Robert-Asselin filter parameter and the viscous and diffusive coefficients.
-
+numerical diffusion in both the time and space discretisations of the 
+advective term: they are both set explicitly by the user through the Robert-Asselin 
+filter parameter and the viscous and diffusive coefficients.
+
+\gmcomment{
 %gm - reflexion about leapfrog: ongoing work with Matthieu Leclair
 % to be updated latter with addition of new time stepping strategy
-\amtcomment{
 \colorbox{yellow}{Note}: 
-1- There is no reason why one should apply a same value of $\gamma$ on both momentum and tracer equations. In climate applications, one could found useful to use a lower value on tracer (quantity that one wants to conserve) than on the dynamics. We never explore this possibility.
-The Robert-Asselin time filter slightly departs from a simple second order time diffusive operator computed with a forward time stepping due to the presence of $x_f^{t-\Delta t}$ in the right hand side of  \ref{Eq_DOM_nxt_asselin}. The original willing of Robert1966 and Asselin1972 was to design a time filter that allow much larger parameter than 0.5.   is due to computer saving consideration. In the original asselin filter, $x^{t-\Delta t}$ is used instead:
+The Robert-Asselin time filter slightly departs from a simple second order time 
+diffusive operator computed with a forward time stepping due to the presence of 
+$x_f^{t-\Delta t}$ in the right hand side of  \ref{Eq_DOM_nxt_asselin}. The original 
+willing of Robert1966 and Asselin1972 was to design a time filter that allow much 
+larger parameter than 0.5.   is due to computer saving consideration. In the original 
+asselin filter, $x^{t-\Delta t}$ is used instead:
  \begin{equation} \label{Eq_DOM_nxt_asselin_true}
 x_f^t  = x^t + \gamma \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right]
 \end{equation} 
-Applying a "true" Asselin time filter is nothing more than adding a harmonic diffusive operator in time. Indeed, equations \ref{Eq_DOM_nxt} and \ref{Eq_DOM_nxt_asselin_true} can be rewritten together as:
+Applying a "true" Asselin time filter is nothing more than adding a harmonic 
+diffusive operator in time. Indeed, equations \ref{Eq_DOM_nxt} and 
+\ref{Eq_DOM_nxt_asselin_true} can be rewritten together as:
 \begin{equation} \label{Eq_DOM_nxt2}
 \begin{split}
   \frac{ x^{t+\Delta t} - x^{t-\Delta t} } { 2 \,\Delta t } 
   &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} + \frac{ x_f^t  - x^t }{2 \,\Delta t} \\
-  &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} + \gamma\ \frac{  \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] }{2 \,\Delta t}  \\
+  &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} 
+    + \gamma\ \frac{  \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] }{2 \,\Delta t}  \\
   &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} 
   + 2 \Delta t \ \gamma \ \frac{1}{{2 \Delta t}^2}  
@@ -621 +950 @@
 \end{equation} 
 
-Equations  \ref{Eq_DOM_nxt2} and \ref{Eq_DOM_nxt3} suggest several remarks. First the Asselin filter is definitively a second order time diffusive operator which is evaluated at centered time step. The magnitude of this diffusion is proportional to the time step (with $\gamma$ usually taken between $10^{-1}$ to $10^{-3}$) . Second, this term has to be taken into account in all budgets of the equations (mass, heat content, salt content, kinetic energy...). Nevertheless,we stress here that it is small and does not create systematic biases. Indeed let us evaluate how it contributes to the time evolution of $x$ between $t_o$ and $t_1$:
+Equations  \ref{Eq_DOM_nxt2} and \ref{Eq_DOM_nxt3} suggest several remarks. 
+First the Asselin filter is definitively a second order time diffusive operator which is 
+evaluated at centered time step. The magnitude of this diffusion is proportional to 
+the time step (with $\gamma$ usually taken between $10^{-1}$ to $10^{-3}$). 
+Second, this term has to be taken into account in all budgets of the equations 
+(mass, heat content, salt content, kinetic energy...). Nevertheless,we stress here 
+that it is small and does not create systematic biases. Indeed let us evaluate how 
+it contributes to the time evolution of $x$ between $t_o$ and $t_1$:
 \begin{equation} \label{Eq_DOM_nxt4}
 \begin{split}
@@ -640 +976 @@
 \label{DOM_nxt_forward_imp}
 
-The leapfrog differencing is unsuitable for the representation of diffusive 
-and damping processes. For $D$, a horizontal diffusive terms and/or the 
-restoring terms to a tracer climatology (when they are present, see 
-\S~\ref{TRA_dmp}), a forward time differencing scheme is used :
+The leapfrog differencing scheme is unsuitable for the representation of 
+diffusive and damping processes. For a tendancy $D_x$, representing a 
+diffusive term or a restoring term to a tracer climatology 
+(when present, see \S~\ref{TRA_dmp}), a forward time differencing scheme
+ is used :
 \begin{equation} \label{Eq_DOM_nxt_euler}
-   x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t-\Delta t}
+   x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ {D_x}^{t-\Delta t}
 \end{equation} 
 
 This is diffusive in time and conditionally stable. For example, the 
-condition of stability for a second and fourth order horizontal diffusions are \citep{Griffies2004}:
+conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}:
 \begin{equation} \label{Eq_DOM_nxt_euler_stability}
 A^h < \left\{
@@ -658 +995 @@
 \right.
 \end{equation}
-where $e$ is the smallest grid size in the two horizontal direction and $A^h$ 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 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 
 used, but usually the numerical stability condition implies a strong 
-constraint on the time step. Two solutions are available in OPA to overcome 
+constraint on the time step. Two solutions are available in \NEMO to overcome 
 the stability constraint: $(a)$ a forward time differencing scheme using a 
-time splitting technique (\np{ln\_zdfexp}=T) or $(b)$ a backward (or implicit) 
-time differencing scheme by \np{ln\_zdfexp}=F). In $(a)$, the master 
+time splitting technique (\np{ln\_zdfexp}=.true.) or $(b)$ a backward (or implicit) 
+time differencing scheme by \np{ln\_zdfexp}=.false.). In $(a)$, the master 
 time step $\Delta $t is cut into $N$ fractional time steps so that the 
 stability criterion is reduced by a factor of $N$. The computation is done as 
@@ -678 +1015 @@
 \end{split}
 \end{equation}
-with DF a vertical diffusion term. The number of fractional time steps, $N$, is given by setting \np{n\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally stable but diffusive. It can be written as follows:
+with DF a vertical diffusion term. The number of fractional time steps, $N$, is given 
+by setting \np{n\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally 
+stable but diffusive. It can be written as follows:
 \begin{equation} \label{Eq_DOM_nxt_imp}
    x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t+\Delta t}
@@ -692 +1031 @@
 \right]
 \end{equation}
-where RHS is the right hand side of the equation except 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)
@@ -703 +1042 @@
 \end{align*}
 
-\eqref{Eq_DOM_nxt_imp_mat} is a linear system of equations. All the elements of the corresponding matrix vanish except those on the diagonals. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the 
-two extra-diagonal terms, therefore a special adaptation of the Gauss elimination procedure is used to find the solution (see for example \citet{Richtmyer1967}).
+\eqref{Eq_DOM_nxt_imp_mat} is a linear system of equations which associated 
+matrix is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal 
+term is greater than the sum of the two extra-diagonal terms, therefore a special 
+adaptation of the Gauss elimination procedure is used to find the solution 
+(see for example \citet{Richtmyer1967}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -714 +1056 @@
 \namdisplay{namrun}         
 %--------------------------------------------------------------------------------------------------------------
-The first time step of this three level scheme when starting from initial conditions is a forward step (Euler time integration): $x^1 = x^0 + \Delta t \ \text{RHS}^0$.
+
+The first time step of this three level scheme when starting from initial conditions 
+is a forward step (Euler time integration): 
+\begin{equation} \label{Eq_DOM_euler}
+   x^1 = x^0 + \Delta t \ \text{RHS}^0
+\end{equation}
 
 It is also possible to restart from a previous computation, by using a 
@@ -720 +1067 @@
 restartability of the code: the user should obtain the same results to 
 machine precision either by running the model for $2N$ time steps in one go, 
-or by performing two consecutive experiments of $N$ steps with a restart. This 
-requires saving two time levels and many auxiliary data in the restart files 
-in double precision. 
-
+or by performing two consecutive experiments of $N$ steps with a restart. 
+This requires saving two time levels and many auxiliary data in the restart 
+files in machine precision. 
+
+Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure 
+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).
+%%%
+\gmcomment{add here how to force the restart to contain only one time step for operational purposes}
+%%%
+
+\gmcomment{       % add a subsection here  
 
 %-------------------------------------------------------------------------------------------------------------
@@ -735 +1091 @@
  \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj}
 
-
-
+}        %% end add
+