Index: /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 11334)
+++ /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 11335)
@@ -2991,4 +2991,15 @@
 }
 
+@article{         white.hoskins.ea_QJRMS05,
+  title         = "Consistent approximate models of the global atmosphere: shallow, deep,
+                  hydrostatic, quasi-hydrostatic and non-hydrostatic",
+  pages         = "2081--2107",
+  journal       = "Quarterly Journal of the Royal Meteorological Society",
+  volume        = "131",
+  author        = "A. A. White and B. J. Hoskins and I. Roulstone and A. Staniforth",
+  year          = "2005",
+  doi           = "10.1256/qj.04.49"
+}
+
 @article{         white.adcroft.ea_JCP09,
   title         = "High-order regridding-remapping schemes for continuous
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex	(revision 11334)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/annex_A.tex	(revision 11335)
@@ -29,38 +29,85 @@
 \begin{equation}
   \label{apdx:A_s_slope}
-  \sigma_1 =\frac{1}{e_1 }\;\left. {\frac{\partial z}{\partial i}} \right|_s
+  \sigma_1 =\frac{1}{e_1 } \; \left. {\frac{\partial z}{\partial i}} \right|_s
   \quad \text{and} \quad
-  \sigma_2 =\frac{1}{e_2 }\;\left. {\frac{\partial z}{\partial j}} \right|_s
-\end{equation}
-
-The chain rule to establish the model equations in the curvilinear $s-$coordinate system is:
+  \sigma_2 =\frac{1}{e_2 } \; \left. {\frac{\partial z}{\partial j}} \right|_s .
+\end{equation}
+
+The model fields (e.g. pressure $p$) can be viewed as functions of $(i,j,z,t)$ (e.g. $p(i,j,z,t)$) or as
+functions of $(i,j,s,t)$ (e.g. $p(i,j,s,t)$). The symbol $\bullet$ will be used to represent any one of 
+these fields.  Any ``infinitesimal'' change in $\bullet$ can be written in two forms: 
+\begin{equation}
+  \label{apdx:A_s_infin_changes}
+  \begin{aligned}
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,s,t} 
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,s,t} 
+                + \delta s \left. \frac{ \partial \bullet }{\partial s} \right|_{i,j,t} 
+                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,s} , \\
+    & \delta \bullet =  \delta i \left. \frac{ \partial \bullet }{\partial i} \right|_{j,z,t} 
+                + \delta j \left. \frac{ \partial \bullet }{\partial i} \right|_{i,z,t} 
+                + \delta z \left. \frac{ \partial \bullet }{\partial z} \right|_{i,j,t} 
+                + \delta t \left. \frac{ \partial \bullet }{\partial t} \right|_{i,j,z} .
+  \end{aligned}
+\end{equation}
+Using the first form and considering a change $\delta i$ with $j, z$ and $t$ held constant, shows that
+\begin{equation}
+  \label{apdx:A_s_chain_rule}
+      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,z,t}  =
+      \left. {\frac{\partial \bullet }{\partial i}} \right|_{j,s,t}
+    + \left. {\frac{\partial s       }{\partial i}} \right|_{j,z,t} \; 
+      \left. {\frac{\partial \bullet }{\partial s}} \right|_{i,j,t} .    
+\end{equation}
+The term $\left. \partial s / \partial i \right|_{j,z,t}$ can be related to the slope of constant $s$ surfaces, 
+(\autoref{apdx:A_s_slope}), by applying the second of (\autoref{apdx:A_s_infin_changes}) with $\bullet$ set to 
+$s$ and $j, t$ held constant
+\begin{equation}
+\label{apdx:a_delta_s}
+\delta s|_{j,t} = 
+         \delta i \left. \frac{ \partial s }{\partial i} \right|_{j,z,t} 
+       + \delta z \left. \frac{ \partial s }{\partial z} \right|_{i,j,t} .
+\end{equation}
+Choosing to look at a direction in the $(i,z)$ plane in which $\delta s = 0$ and using
+(\autoref{apdx:A_s_slope}) we obtain 
+\begin{equation}
+\left. \frac{ \partial s }{\partial i} \right|_{j,z,t} =  
+         -  \left. \frac{ \partial z }{\partial i} \right|_{j,s,t} \;
+            \left. \frac{ \partial s }{\partial z} \right|_{i,j,t}
+	 = - \frac{e_1 }{e_3 }\sigma_1  .
+\label{apdx:a_ds_di_z}
+\end{equation}
+Another identity, similar in form to (\autoref{apdx:a_ds_di_z}), can be derived
+by choosing $\bullet$ to be $s$ and using the second form of (\autoref{apdx:A_s_infin_changes}) to consider
+changes in which $i , j$ and $s$ are constant. This shows that
+\begin{equation}
+\label{apdx:A_w_in_s}
+w_s = \left. \frac{ \partial z }{\partial t} \right|_{i,j,s} =  
+- \left. \frac{ \partial z }{\partial s} \right|_{i,j,t}
+  \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} 
+  = - e_3 \left. \frac{ \partial s }{\partial t} \right|_{i,j,z} . 
+\end{equation}
+
+In what follows, for brevity, indication of the constancy of the $i, j$ and $t$ indices is 
+usually omitted. Using the arguments outlined above one can show that the chain rules needed to establish 
+the model equations in the curvilinear $s-$coordinate system are:
 \begin{equation}
   \label{apdx:A_s_chain_rule}
   \begin{aligned}
     &\left. {\frac{\partial \bullet }{\partial t}} \right|_z  =
-    \left. {\frac{\partial \bullet }{\partial t}} \right|_s
-	 -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial t} \\
+    \left. {\frac{\partial \bullet }{\partial t}} \right|_s 
+	 + \frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial t} , \\
     &\left. {\frac{\partial \bullet }{\partial i}} \right|_z  =
     \left. {\frac{\partial \bullet }{\partial i}} \right|_s
-    -\frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial i}=
-    \left. {\frac{\partial \bullet }{\partial i}} \right|_s
-    -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} \\
+    +\frac{\partial \bullet }{\partial s}\; \frac{\partial s}{\partial i}=
+    \left. {\frac{\partial \bullet }{\partial i}} \right|_s 
+    -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial \bullet }{\partial s} , \\
     &\left. {\frac{\partial \bullet }{\partial j}} \right|_z  =
-    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
-    - \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=
-    \left. {\frac{\partial \bullet }{\partial j}} \right|_s
-    - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} \\
-    &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s}
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
+    + \frac{\partial \bullet }{\partial s}\;\frac{\partial s}{\partial j}=
+    \left. {\frac{\partial \bullet }{\partial j}} \right|_s 
+    - \frac{e_2 }{e_3 }\sigma_2 \frac{\partial \bullet }{\partial s} , \\
+    &\;\frac{\partial \bullet }{\partial z}  \;\; = \frac{1}{e_3 }\frac{\partial \bullet }{\partial s} .
   \end{aligned}
 \end{equation}
 
-In particular applying the time derivative chain rule to $z$ provides the expression for $w_s$,
-the vertical velocity of the $s-$surfaces referenced to a fix z-coordinate:
-\begin{equation}
-  \label{apdx:A_w_in_s}
-  w_s   =  \left.   \frac{\partial z }{\partial t}   \right|_s
-  = \frac{\partial z}{\partial s} \; \frac{\partial s}{\partial t}
-  = e_3 \, \frac{\partial s}{\partial t}
-\end{equation}
 
 % ================================================================
@@ -131,43 +178,42 @@
 \end{subequations}
 
-Here, $w$ is the vertical velocity relative to the $z-$coordinate system.
-Introducing the dia-surface velocity component,
-$\omega $, defined as the volume flux across the moving $s$-surfaces per unit horizontal area:
+Here, $w$ is the vertical velocity relative to the $z-$coordinate system. 
+Using the first form of (\autoref{apdx:A_s_infin_changes}) 
+and the definitions (\autoref{apdx:A_s_slope}) and (\autoref{apdx:A_w_in_s}) for $\sigma_1$, $\sigma_2$ and  $w_s$,
+one can show that the vertical velocity, $w_p$ of a point
+moving with the horizontal velocity of the fluid along an $s$ surface is given by 
+\begin{equation}
+\label{apdx:A_w_p}
+\begin{split}
+w_p  = & \left. \frac{ \partial z }{\partial t} \right|_s
+     + \frac{u}{e_1} \left. \frac{ \partial z }{\partial i} \right|_s 
+     + \frac{v}{e_2} \left. \frac{ \partial z }{\partial j} \right|_s \\
+     = & w_s + u \sigma_1 + v \sigma_2 .
+\end{split}     
+\end{equation}
+ The vertical velocity across this surface is denoted by
 \begin{equation}
   \label{apdx:A_w_s}
-  \omega  = w - w_s - \sigma_1 \,u - \sigma_2 \,v    \\
-\end{equation}
-with $w_s$ given by \autoref{apdx:A_w_in_s},
-we obtain the expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
-\begin{subequations}
-  \begin{align*}
-    {
-    \begin{array}{*{20}l}
-      \nabla \cdot {\mathrm {\mathbf U}}
-      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
+  \omega  = w - w_p = w - ( w_s + \sigma_1 \,u + \sigma_2 \,v )  . 
+\end{equation}
+Hence 
+\begin{equation}
+\frac{1}{e_3 } \frac{\partial}{\partial s}   \left[  w -  u\;\sigma_1  - v\;\sigma_2  \right] = 
+\frac{1}{e_3 } \frac{\partial}{\partial s} \left[  \omega + w_s \right] = 
+   \frac{1}{e_3 } \left[ \frac{\partial \omega}{\partial s} 
+ + \left. \frac{ \partial }{\partial t} \right|_s \frac{\partial z}{\partial s} \right] = 
+   \frac{1}{e_3 } \frac{\partial \omega}{\partial s} + \frac{1}{e_3 } \left. \frac{ \partial e_3}{\partial t} . \right|_s 
+\end{equation}
+
+Using (\autoref{apdx:A_w_s}) in our expression for $\nabla \cdot {\mathrm {\mathbf U}}$ we obtain 
+our final expression for the divergence of the velocity in the curvilinear $s-$coordinate system:
+\begin{equation}
+      \nabla \cdot {\mathrm {\mathbf U}} =
+         \frac{1}{e_1 \,e_2 \,e_3 }    \left[
 		  \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
         +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
         + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
-        + \frac{1}{e_3 } \frac{\partial w_s       }{\partial s} \\ \\
-      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
-		  \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
-        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
-        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
-        + \frac{1}{e_3 } \frac{\partial}{\partial s}  \left(  e_3 \; \frac{\partial s}{\partial t}   \right) \\ \\
-      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
-		  \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
-        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
-        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
-        + \frac{\partial}{\partial s} \frac{\partial s}{\partial t}
-        + \frac{1}{e_3 } \frac{\partial s}{\partial t} \frac{\partial e_3}{\partial s} \\ \\
-      &= \frac{1}{e_1 \,e_2 \,e_3 }    \left[
-		  \left.  \frac{\partial (e_2 \,e_3 \,u)}{\partial i} \right|_s
-        +\left.  \frac{\partial (e_1 \,e_3 \,v)}{\partial j} \right|_s        \right]
-        + \frac{1}{e_3 } \frac{\partial \omega }{\partial s}
-        + \frac{1}{e_3 } \frac{\partial e_3}{\partial t}
-    \end{array}
-        }
-  \end{align*}
-\end{subequations}
+        + \frac{1}{e_3 } \left. \frac{\partial e_3}{\partial t} \right|_s .
+\end{equation}
 
 As a result, the continuity equation \autoref{eq:PE_continuity} in the $s-$coordinates is:
@@ -178,7 +224,7 @@
     {\left. {\frac{\partial (e_2 \,e_3 \,u)}{\partial i}} \right|_s
       +  \left. {\frac{\partial (e_1 \,e_3 \,v)}{\partial j}} \right|_s } \right]
-  +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0
-\end{equation}
-A additional term has appeared that take into account
+  +\frac{1}{e_3 }\frac{\partial \omega }{\partial s} = 0 .
+\end{equation}
+An additional term has appeared that takes into account
 the contribution of the time variation of the vertical coordinate to the volume budget.
 
@@ -210,6 +256,5 @@
         + w \;\frac{\partial u}{\partial z} \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
-        - \left. \zeta \right|_z v
-        +  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z
+        -  \frac{1}{e_1 \,e_2 }\left[ { \left.{ \frac{\partial (e_2 \,v)}{\partial i} }\right|_z
         -\left.{ \frac{\partial (e_1 \,u)}{\partial j} }\right|_z } \right] \; v
         +  \frac{1}{2e_1} \left.{ \frac{\partial (u^2+v^2)}{\partial i} } \right|_z
@@ -230,13 +275,13 @@
         } \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
-        + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\
       &\qquad \qquad \qquad \quad
         + \frac{w}{e_3 } \;\frac{\partial u}{\partial s}
-        - \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s}
+        + \left[   {\frac{\sigma_1 }{e_3 }\frac{\partial v}{\partial s}
         - \frac{\sigma_2 }{e_3 }\frac{\partial u}{\partial s}}   \right]\;v
         - \frac{\sigma_1 }{2e_3 }\frac{\partial (u^2+v^2)}{\partial s} \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
-        + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s \\
       &\qquad \qquad \qquad \quad
@@ -245,15 +290,15 @@
         - \sigma_1 u\frac{\partial u}{\partial s} - \sigma_1 v\frac{\partial v}{\partial s}} \right] \\ \\
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
-        + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
         + \frac{1}{e_3} \left[  w - \sigma_2 v - \sigma_1 u  \right]
-        \; \frac{\partial u}{\partial s}   \\
+        \; \frac{\partial u}{\partial s} .  \\
         %
-      \intertext{Introducing $\omega$, the dia-a-surface velocity given by (\autoref{apdx:A_w_s}) }
+      \intertext{Introducing $\omega$, the dia-s-surface velocity given by (\autoref{apdx:A_w_s}) }
       %
       &= \left. {\frac{\partial u }{\partial t}} \right|_z
-        + \left. \zeta \right|_s \;v
+        - \left. \zeta \right|_s \;v
         + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
-        + \frac{1}{e_3 } \left( \omega - w_s \right) \frac{\partial u}{\partial s}   \\
+        + \frac{1}{e_3 } \left( \omega + w_s \right) \frac{\partial u}{\partial s}   \\
     \end{array}
     }
@@ -266,11 +311,11 @@
   {
     \begin{array}{*{20}l}
-      w_s  \;\frac{\partial u}{\partial s}
-      = \frac{\partial s}{\partial t} \;  \frac{\partial u }{\partial s}
-      = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \quad ,
+      \frac{w_s}{e_3}  \;\frac{\partial u}{\partial s}
+      = - \left. \frac{\partial s}{\partial t} \right|_z \;  \frac{\partial u }{\partial s}
+      = \left. {\frac{\partial u }{\partial t}} \right|_s  - \left. {\frac{\partial u }{\partial t}} \right|_z \ .
     \end{array}
   }
 \]
-leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,
+This leads to the $s-$coordinate formulation of the total $z-$coordinate time derivative,
 \ie the total $s-$coordinate time derivative :
 \begin{align}
@@ -278,7 +323,7 @@
   \left. \frac{D u}{D t} \right|_s
   = \left. {\frac{\partial u }{\partial t}} \right|_s
-  + \left. \zeta \right|_s \;v
+  - \left. \zeta \right|_s \;v
   + \frac{1}{2\,e_1}\left. {\frac{\partial (u^2+v^2)}{\partial i}} \right|_s
-  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s}
+  + \frac{1}{e_3 } \omega \;\frac{\partial u}{\partial s} . 
 \end{align}
 Therefore, the vector invariant form of the total time derivative has exactly the same mathematical form in
@@ -306,5 +351,5 @@
                                          + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
                                       &&- \frac{v}{e_1 e_2 }\left(    v	\;\frac{\partial e_2 }{\partial i}
-                                         -u	\;\frac{\partial e_1 }{\partial j} 	\right) \\
+                                         -u	\;\frac{\partial e_1 }{\partial j} 	\right) . \\
   \end{array}
   }
@@ -328,5 +373,5 @@
        -  e_2 u \;\frac{\partial e_3 }{\partial i}
        -  e_1 v \;\frac{\partial e_3 }{\partial j}   \right)
-       -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
+       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right] \\ \\
     && - \frac{v}{e_1 e_2 }\left( 	v  \;\frac{\partial e_2 }{\partial i}
        -u  \;\frac{\partial e_1 }{\partial j} 	\right) \\ \\
@@ -337,7 +382,7 @@
     && - \,u \left[  \frac{1}{e_1 e_2 e_3} \left(   \frac{\partial(e_2 e_3 \, u)}{\partial i}
        + \frac{\partial(e_1 e_3 \, v)}{\partial j}  \right)
-       -\frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right]
+       + \frac{1}{e_3}        \frac{\partial \omega}{\partial s}                       \right]
        - \frac{v}{e_1 e_2 }\left( 	v   \;\frac{\partial e_2 }{\partial i}
-       -u   \;\frac{\partial e_1 }{\partial j} 	\right)                  \\
+       -u   \;\frac{\partial e_1 }{\partial j} 	\right)     .             \\
      %
     \intertext {Introducing a more compact form for the divergence of the momentum fluxes,
@@ -361,5 +406,5 @@
   + \left.  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)    \right|_s
   - \frac{v}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
-    -u  \;\frac{\partial e_1 }{\partial j}            \right)
+    -u  \;\frac{\partial e_1 }{\partial j}            \right).
 \end{flalign}
 which is the total time derivative expressed in the curvilinear $s-$coordinate system.
@@ -377,9 +422,9 @@
     & =-\frac{1}{\rho_o e_1 }\left[ {\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{e_1 }{e_3 }\sigma_1 \frac{\partial p}{\partial s}} \right] \\
     & =-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s +\frac{\sigma_1 }{\rho_o \,e_3 }\left( {-g\;\rho \;e_3 } \right) \\
-    &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1
+    &=-\frac{1}{\rho_o \,e_1 }\left. {\frac{\partial p}{\partial i}} \right|_s -\frac{g\;\rho }{\rho_o }\sigma_1 .
   \end{split}
 \]
 Applying similar manipulation to the second component and
-replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it comes:
+replacing $\sigma_1$ and $\sigma_2$ by their expression \autoref{apdx:A_s_slope}, it becomes:
 \begin{equation}
   \label{apdx:A_grad_p_1}
@@ -391,5 +436,5 @@
     -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z
     &=-\frac{1}{\rho_o \,e_2 } \left(    \left.               {\frac{\partial p}{\partial j}} \right|_s
-      + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) \\
+      + g\;\rho \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right) . \\
   \end{split}
 \end{equation}
@@ -405,6 +450,6 @@
 \[
   \begin{split}
-    p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \left(  \rho_o \, d + 1 \right) \; e_3 \; dk   \\
-    &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + g \, \int_z^\eta e_3 \; dk
+    p &= g\; \int_z^\eta \rho \; e_3 \; dk = g\; \int_z^\eta \rho_o \left( d + 1 \right) \; e_3 \; dk   \\
+    &= g \, \rho_o \; \int_z^\eta d \; e_3 \; dk + \rho_o g \, \int_z^\eta e_3 \; dk .
   \end{split}
 \]
@@ -412,13 +457,13 @@
 \begin{equation}
   \label{apdx:A_pressure}
-  p = \rho_o \; p_h' + g \, ( z + \eta )
+  p = \rho_o \; p_h' + \rho_o \, g \, ( \eta - z )
 \end{equation}
 and the hydrostatic pressure balance expressed in terms of $p_h'$ and $d$ is:
 \[
-  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3
+  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 .
 \]
 
 Substituing \autoref{apdx:A_pressure} in \autoref{apdx:A_grad_p_1} and
-using the definition of the density anomaly it comes the expression in two parts:
+using the definition of the density anomaly it becomes an expression in two parts:
 \begin{equation}
   \label{apdx:A_grad_p_2}
@@ -426,9 +471,9 @@
     -\frac{1}{\rho_o \, e_1 } \left. {\frac{\partial p}{\partial i}} \right|_z
     &=-\frac{1}{e_1 } \left(     \left.              {\frac{\partial p_h'}{\partial i}} \right|_s
-      + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} \\
+      + g\; d  \;\left. {\frac{\partial z}{\partial i}} \right|_s    \right)  - \frac{g}{e_1 } \frac{\partial \eta}{\partial i} ,  \\
              %
     -\frac{1}{\rho_o \, e_2 }\left. {\frac{\partial p}{\partial j}} \right|_z
     &=-\frac{1}{e_2 } \left(    \left.               {\frac{\partial p_h'}{\partial j}} \right|_s
-      + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j}\\
+      + g\; d \;\left. {\frac{\partial z}{\partial j}} \right|_s   \right)  - \frac{g}{e_2 } \frac{\partial \eta}{\partial j} . \\
   \end{split}
 \end{equation}
@@ -463,5 +508,5 @@
     -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right)
     -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i}
-    +   D_u^{\vect{U}}  +   F_u^{\vect{U}}
+    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} ,
   \end{multline}
   \begin{multline}
@@ -473,5 +518,5 @@
     -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right)
     -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j}
-    +  D_v^{\vect{U}}  +   F_v^{\vect{U}}
+    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} .
   \end{multline}
 \end{subequations}
@@ -483,10 +528,10 @@
     \label{apdx:A_PE_dyn_flux_u}
     \frac{1}{e_3} \frac{\partial \left(  e_3\,u  \right) }{\partial t} =
-    \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)
+    - \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,u}   \right)
     +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
           -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,v     \\
     -   \frac{1}{e_1 } \left(    \frac{\partial p_h'}{\partial i} + g\; d  \; \frac{\partial z}{\partial i}    \right)
     -   \frac{g}{e_1 } \frac{\partial \eta}{\partial i}
-    +   D_u^{\vect{U}}  +   F_u^{\vect{U}}
+    +   D_u^{\vect{U}}  +   F_u^{\vect{U}} ,
   \end{multline}
   \begin{multline}
@@ -494,9 +539,9 @@
     \frac{1}{e_3}\frac{\partial \left(  e_3\,v  \right) }{\partial t}=
     -  \nabla \cdot \left(   {{\mathrm {\mathbf U}}\,v}   \right)
-    +   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
+    -   \left\{ {f + \frac{1}{e_1 e_2 }\left(    v  \;\frac{\partial e_2 }{\partial i}
           -u  \;\frac{\partial e_1 }{\partial j}            \right)} \right\} \,u     \\
     -   \frac{1}{e_2 } \left(    \frac{\partial p_h'}{\partial j} + g\; d  \; \frac{\partial z}{\partial j}    \right)
     -   \frac{g}{e_2 } \frac{\partial \eta}{\partial j}
-    +  D_v^{\vect{U}}  +   F_v^{\vect{U}}
+    +  D_v^{\vect{U}}  +   F_v^{\vect{U}} . 
   \end{multline}
 \end{subequations}
@@ -505,5 +550,5 @@
 \begin{equation}
   \label{apdx:A_dyn_zph}
-  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3
+  \frac{\partial p_h'}{\partial k} = - d \, g \, e_3 .
 \end{equation}
 
@@ -531,9 +576,9 @@
   \left[           \frac{\partial }{\partial i} \left( {e_2 \,e_3 \;Tu} \right)
     +   \frac{\partial }{\partial j} \left( {e_1 \,e_3 \;Tv} \right)               \right]       \\
-  +  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right)
+  -  \frac{1}{e_3}  \frac{\partial }{\partial k} \left(                   Tw  \right)
   +  D^{T} +F^{T}
 \end{multline}
 
-The expression for the advection term is a straight consequence of (A.4),
+The expression for the advection term is a straight consequence of (\autoref{apdx:A_sco_Continuity}),
 the expression of the 3D divergence in the $s-$coordinates established above. 
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex	(revision 11334)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/annex_B.tex	(revision 11335)
@@ -53,9 +53,14 @@
   {
   \begin{array}{*{20}l}
-    D^T=& \frac{1}{e_1\,e_2\,e_3 }\;\left[ {\ \ \ \ e_2\,e_3\,A^{lT} \;\left.
-          {\frac{\partial }{\partial i}\left( {\frac{1}{e_1}\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_1 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.  \\
-        &\qquad \quad \ \ \ +e_1\,e_3\,A^{lT} \;\left. {\frac{\partial }{\partial j}\left( {\frac{1}{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s -\frac{\sigma_2 }{e_3 }\;\frac{\partial T}{\partial s}} \right)} \right|_s \\
-        &\qquad \quad \ \ \ + e_1\,e_2\,A^{lT} \;\frac{\partial }{\partial s}\left( {-\frac{\sigma_1 }{e_1 }\;\left. {\frac{\partial T}{\partial i}} \right|_s -\frac{\sigma_2 }{e_2 }\;\left. {\frac{\partial T}{\partial j}} \right|_s } \right.
-          \left. {\left. {+\left( {\varepsilon +\sigma_1^2+\sigma_2 ^2} \right)\;\frac{1}{e_3 }\;\frac{\partial T}{\partial s}} \right)\;\;} \right]
+    D^T= \frac{1}{e_1\,e_2\,e_3 } & \left\{ \quad \quad \frac{\partial }{\partial i}  \left. \left[  e_2\,e_3 \, A^{lT}
+                               \left( \  \frac{1}{e_1}\; \left. \frac{\partial T}{\partial i} \right|_s 
+                                       -\frac{\sigma_1 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
+        &  \quad \  +   \            \left.   \frac{\partial }{\partial j}  \left. \left[  e_1\,e_3 \, A^{lT}
+                               \left( \ \frac{1}{e_2 }\; \left. \frac{\partial T}{\partial j} \right|_s 
+                                       -\frac{\sigma_2 }{e_3 } \; \frac{\partial T}{\partial s} \right) \right]  \right|_s  \right. \\
+        &  \quad \  +   \           \left.  e_1\,e_2\, \frac{\partial }{\partial s}  \left[ A^{lT} \; \left( 
+	                  -\frac{\sigma_1 }{e_1 } \; \left. \frac{\partial T}{\partial i} \right|_s 
+	                  -\frac{\sigma_2 }{e_2 } \; \left. \frac{\partial T}{\partial j} \right|_s 
+                          +\left( \varepsilon +\sigma_1^2+\sigma_2 ^2 \right) \; \frac{1}{e_3 } \; \frac{\partial T}{\partial s} \right) \; \right] \;  \right\} .
   \end{array}
           }
@@ -90,7 +95,7 @@
   {
   \begin{array}{*{20}l}
-    \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, it becomes:}
+    \intertext{Noting that $\frac{1}{e_1} \left. \frac{\partial e_3 }{\partial i} \right|_s = \frac{\partial \sigma_1 }{\partial s}$, this becomes:}
     %
-    & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
+    D^T & =\frac{1}{e_1\,e_2\,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2\,e_3 }{e_1}\,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)} \right|_s \left. -\, {e_3 \frac{\partial }{\partial i}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
     & \qquad \qquad \quad -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s -e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\
     & \qquad \qquad \quad\shoveright{ \left. { +e_1 \,\sigma_1 \frac{\partial }{\partial s}\left( {\frac{e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }\\
@@ -99,5 +104,5 @@
     & \qquad \qquad \quad \left. {+\frac{e_2 \,\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s} \;\frac{\partial e_3 }{\partial i}}  \right|_s -e_2 A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s \\
     & \qquad \qquad \quad-e_2 \,\sigma_1 \frac{\partial}{\partial s}\left( {A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}} \right) \\
-    & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]}
+    & \qquad \qquad \quad\shoveright{ \left. {-\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s} \left( {\frac{\sigma_1 }{e_3}A^{lT}\;\frac{\partial T}{\partial s}} \right) + \frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 }{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} .
   \end{array}
       } \\
@@ -105,10 +110,10 @@
   \begin{array}{*{20}l}
     %
-    \intertext{Using the same remark as just above, it becomes:}
+    \intertext{Using the same remark as just above, $D^T$ becomes:}
     %
-    &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\
+   D^T &= \frac{1}{e_1 \,e_2 \,e_3 } \left[ {\left. {\;\;\;\frac{\partial }{\partial i} \left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right.\;\;\; \\
     & \qquad \qquad \quad+\frac{e_1 \,e_2 \,\sigma_1 }{e_3 }A^{lT}\;\frac{\partial T}{\partial s}\;\frac{\partial \sigma_1 }{\partial s} - \frac {\sigma_1 }{e_3} A^{lT} \;\frac{\partial \left( {e_1 \,e_2 \,\sigma_1 } \right)}{\partial s}\;\frac{\partial T}{\partial s} \\
     & \qquad \qquad \quad-e_2 \left( {A^{lT}\;\frac{\partial \sigma_1 }{\partial s}\left. {\frac{\partial T}{\partial i}} \right|_s +\frac{\partial }{\partial s}\left( {\sigma_1 A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right)-\frac{\partial \sigma_1 }{\partial s}\;A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s } \right) \\
-    & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] }
+    & \qquad \qquad \quad\shoveright{\left. {+\frac{\partial }{\partial s}\left( {\frac{e_1 \,e_2 \,\sigma_1 ^2}{e_3 }A^{lT}\;\frac{\partial T}{\partial s}+\frac{e_1 \,e_2}{e_3 }A^{vT}\;\frac{\partial T}{\partial s}} \right)\;\;\;} \right] . }
   \end{array}
       } \\
@@ -120,10 +125,10 @@
     the third line reduces to a single vertical derivative, so it becomes:}
   %
-    & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
+    D^T & =\frac{1}{e_1 \,e_2 \,e_3 }\left[ {\left. {\;\;\;\frac{\partial }{\partial i}\left( {\frac{e_2 \,e_3 }{e_1 }A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s -e_2 \,\sigma_1 \,A^{lT}\;\frac{\partial T}{\partial s}} \right)} \right|_s } \right. \\
     & \qquad \qquad \quad \shoveright{ \left. {+\frac{\partial }{\partial s}\left( {-e_2 \,\sigma_1 \,A^{lT}\;\left. {\frac{\partial T}{\partial i}} \right|_s +A^{lT}\frac{e_1 \,e_2 }{e_3 }\;\left( {\varepsilon +\sigma_1 ^2} \right)\frac{\partial T}{\partial s}} \right)\;\;\;} \right]} \\
     %
     \intertext{In other words, the horizontal/vertical Laplacian operator in the ($i$,$s$) plane takes the following form:}
   \end{array}
-  } \\ 
+  } \\
   %
   {\frac{1}{e_1\,e_2\,e_3}}
@@ -228,6 +233,6 @@
 The isopycnal diffusion operator \autoref{apdx:B4},
 \autoref{apdx:B_ldfiso} conserves tracer quantity and dissipates its square.
-The demonstration of the first property is trivial as \autoref{apdx:B4} is the divergence of fluxes.
-Let us demonstrate the second one:
+As \autoref{apdx:B4} is the divergence of a flux, the demonstration of the first property is trivial, providing that the flux normal to the boundary is zero 
+(as it is when $A_h$ is zero at the boundary). Let us demonstrate the second one:
 \[
   \iiint\limits_D T\;\nabla .\left( {\textbf{A}}_{\textbf{I}} \nabla T \right)\,dv
@@ -248,5 +253,5 @@
              j}-a_2 \frac{\partial T}{\partial k}} \right)^2}
              +\varepsilon \left(\frac{\partial T}{\partial k}\right) ^2\right]      \\
-           & \geq 0
+           & \geq 0 . 
   \end{array}
              }
@@ -365,5 +370,5 @@
 the third component of the second vector is obviously zero and thus :
 \[
-  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right)
+  \Delta {\textbf{U}}_h = \nabla _h \left( \chi \right) - \nabla _h \times \left( \zeta \textbf{k} \right) + \frac {1}{e_3 } \frac {\partial }{\partial k} \left( {\frac {1}{e_3 } \frac{\partial {\textbf{ U}}_h }{\partial k}} \right) . 
 \]
 
@@ -381,5 +386,5 @@
   - \nabla _h \times \left( {A^{lm}\;\zeta \;{\textbf{k}}} \right)
   + \frac{1}{e_3 }\frac{\partial }{\partial k}\left( {\frac{A^{vm}\;}{e_3 }
-      \frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}} \right) \\
+      \frac{\partial {\mathrm {\mathbf U}}_h }{\partial k}} \right) , \\
 \end{equation}
 that is, in expanded form:
@@ -388,9 +393,9 @@
   & = \frac{1}{e_1} \frac{\partial \left( {A^{lm}\chi   } \right)}{\partial i}
     -\frac{1}{e_2} \frac{\partial \left( {A^{lm}\zeta } \right)}{\partial j}
-    +\frac{1}{e_3} \frac{\partial}{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial u}{\partial k} \right)      \\
+    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial u}{\partial k} \right)   ,   \\
   D^{\textbf{U}}_v
   & = \frac{1}{e_2 }\frac{\partial \left( {A^{lm}\chi   } \right)}{\partial j}
     +\frac{1}{e_1 }\frac{\partial \left( {A^{lm}\zeta } \right)}{\partial i}
-    +\frac{1}{e_3} \frac{\partial}{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial v}{\partial k} \right)  
+    +\frac{1}{e_3} \frac{\partial }{\partial k} \left( \frac{A^{vm}}{e_3} \frac{\partial v}{\partial k} \right) .
 \end{align*}
 
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11334)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex	(revision 11335)
@@ -32,6 +32,8 @@
 \begin{enumerate}
 \item
-  \textit{spherical earth approximation}: the geopotential surfaces are assumed to be spheres so that
-  gravity (local vertical) is parallel to the earth's radius
+  \textit{spherical Earth approximation}: the geopotential surfaces are assumed to be oblate spheriods
+  that follow the Earth's bulge; these spheroids are approximated by spheres with
+  gravity locally vertical (parallel to the Earth's radius) and independent of latitude 
+  \citep[][section 2]{white.hoskins.ea_QJRMS05}.   
 \item
   \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius
@@ -63,8 +65,14 @@
     \nabla \cdot \vect U = 0
   \end{equation}
+ \item 
+  \textit{Neglect of additional Coriolis terms}: the Coriolis terms that vary with the cosine of latitude are neglected. 
+  These terms may be non-negligible where the Brunt-Vaisala frequency $N$ is small, either in the deep ocean or
+  in the sub-mesoscale motions of the mixed layer, or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. 
+  They can be consistently included as part of the ocean dynamics \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are 
+  retained in the MIT ocean model.     
 \end{enumerate}
 
 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 $(i,j,k)$ linked to the earth such that
+it is useful to choose an orthogonal set of unit vectors $(i,j,k)$ linked to the Earth such that
 $k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$,
 \ie tangent to the geopotential surfaces.
@@ -107,6 +115,7 @@
 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
+where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface 
+(discretisation can introduce additional artificial ``side-wall'' boundaries). 
+Both $H$ and $\eta$ are referenced to a surface of constant geopotential (\ie a mean sea surface height) on which $z = 0$. 
 (\autoref{fig:ocean_bc}).
 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with
@@ -210,5 +219,5 @@
 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 time step has to be very short when they are present in the model.
 The latter strategy filters out these waves since the rigid lid approximation implies $\eta = 0$,
 \ie the sea surface is the surface $z = 0$.
@@ -217,5 +226,5 @@
 The rigid-lid hypothesis is an obsolescent feature in modern OGCMs.
 It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings.
-Only the free surface formulation is now described in the this document (see the next sub-section).
+Only the free surface formulation is now described in this document (see the next sub-section).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -237,5 +246,5 @@
 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.
+These waves are barotropic (\ie nearly independent of depth) and their phase speed is quite high.
 Their time scale is short with respect to the other processes described by the primitive equations.
 
@@ -266,5 +275,5 @@
 the implicit scheme \citep{dukowicz.smith_JGR94} or
 the addition of a filtering force in the momentum equation \citep{roullet.madec_JGR00}.
-With the present release, \NEMO offers the choice between
+With the present release, \NEMO  offers the choice between
 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or
 a split-explicit scheme strongly inspired the one proposed by \citet{shchepetkin.mcwilliams_OM05}
@@ -338,5 +347,5 @@
 the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$.
 The scalar and vector operators that appear in the primitive equations
-(\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can be written in the tensorial form,
+(\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can then be written in the tensorial form,
 invariant in any orthogonal horizontal curvilinear coordinate system transformation:
 \begin{subequations}
@@ -384,5 +393,5 @@
 \end{gather}
 
-Using the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that
+Using again 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$,
 $NLT$ the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows:
@@ -456,5 +465,5 @@
   &      &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u
             + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\
-  \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it comes:}
+  \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:}
   &      &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v)
 \end{alignat*}
@@ -516,5 +525,5 @@
     % \label{eq:PE_dyn_flux_v}
     \pd[v]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\
-                + \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\
+                - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\
                 - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt)
                 + D_v^{\vect U} + F_v^{\vect U}
@@ -526,5 +535,5 @@
     p_s = \rho \,g \, \eta
   \]
-  with $\eta$ is solution of \autoref{eq:PE_ssh}.
+  and $\eta$ is the solution of \autoref{eq:PE_ssh}.
 
   The vertical velocity and the hydrostatic pressure are diagnosed from the following equations:
@@ -536,15 +545,16 @@
   \]
   where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}.
-\item \textit{tracer equations}:
-  \[
-    %\label{eq:S}
-    \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt]
+
+\item 
+  \textbf{tracer equations}:
+  \begin{equation}
+  \begin{split}
+    \pd[T]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt]
                 - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\
-    %\label{eq:T}
-    \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt]
-                - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S
-    %\label{eq:rho}
-    \rho = \rho \big( T,S,z(k) \big)
-  \]
+    \pd[S]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt]
+                - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\
+    \rho = & \rho \big( T,S,z(k) \big)
+  \end{split}
+  \end{equation}
 \end{itemize}
 
@@ -575,5 +585,5 @@
 follows the isopycnal surfaces, \eg an isopycnic coordinate.
 
-In order to satisfy two or more constrains one can even be tempted to mixed these coordinate systems, as in
+In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, as in
 HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at
 the ocean bottom) \citep{chassignet.smith.ea_JPO03} or
@@ -594,5 +604,5 @@
 This so-called \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact
 an Arbitrary Lagrangian--Eulerian (ALE) coordinate.
-Indeed, choosing an expression for $s$ is an arbitrary choice that determines
+Indeed, one has a great deal of freedom in the choice of expression for $s$. The choice determines
 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and
 which part will be used to move them (Lagrangian part).
@@ -601,13 +611,13 @@
 Its most often used implementation is via an ALE algorithm,
 in which a pure lagrangian step is followed by regridding and remapping steps,
-the later step implicitly embedding the vertical advection
+the latter step implicitly embedding the vertical advection
 \citep{hirt.amsden.ea_JCP74, chassignet.smith.ea_JPO03, white.adcroft.ea_JCP09}.
 Here we follow the \citep{kasahara_MWR74} strategy:
-a regridding step (an update of the vertical coordinate) followed by an eulerian step with
+a regridding step (an update of the vertical coordinate) followed by an Eulerian step with
 an explicit computation of vertical advection relative to the moving s-surfaces.
 
 %\gmcomment{
 %A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient...
-the generalized vertical coordinates used in ocean modelling are not orthogonal,
+The generalized vertical coordinates used in ocean modelling are not orthogonal,
 which contrasts with many other applications in mathematical physics.
 Hence, it is useful to keep in mind the following properties that may seem odd on initial encounter.
@@ -615,5 +625,5 @@
 The horizontal velocity in ocean models measures motions in the horizontal plane,
 perpendicular to the local gravitational field.
-That is, horizontal velocity is mathematically the same regardless the vertical coordinate, be it geopotential,
+That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, be it geopotential,
 isopycnal, pressure, or terrain following.
 The key motivation for maintaining the same horizontal velocity component is that
@@ -660,5 +670,5 @@
 \[
   % \label{eq:PE_sco_w}
-  \omega = w - e_3 \, \pd[s]{t} - \sigma_1 \, u - \sigma_2 \, v
+  \omega = w -  \, \lt. \pd[z]{t} \rt|_s - \sigma_1 \, u - \sigma_2 \, v
 \]
 
@@ -671,5 +681,5 @@
   % \label{eq:PE_sco_u_vector}
     \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} \\
-                - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_1
+                - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1
                 + D_u^{\vect U} + F_u^{\vect U}
   \end{multline*}
@@ -677,5 +687,5 @@
   % \label{eq:PE_sco_v_vector}
     \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j}(u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} \\
-                - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + g \frac{\rho}{\rho_o} \sigma_2
+                - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_2
                 + D_v^{\vect U} + F_v^{\vect U}
   \end{multline*}
@@ -687,5 +697,5 @@
                                        - \frac{1}{e_3} \pd[(\omega \, u)]{k}
                                        - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt)
-                                       + g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U}
+                                       - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U}
   \end{multline*}
   \begin{multline*}
@@ -695,5 +705,5 @@
                                        - \frac{1}{e_3} \pd[(\omega \, v)]{k}
                                        - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt)
-                                       + g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U}
+                                       - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U}
   \end{multline*}
   where the relative vorticity, $\zeta$, the surface pressure gradient,
@@ -750,5 +760,5 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-In that case, the free surface equation is nonlinear, and the variations of volume are fully taken into account.
+In this 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{levier.treguier.ea_rpt07} available on the \NEMO web site.
 
@@ -766,8 +776,13 @@
 The position (\zstar) and vertical discretization (\zstar) are expressed as:
 \[
-  % \label{eq:z-star}
+  % \label{eq:PE_z-star}
   H + \zstar = (H + z)  / r \quad \text{and}  \quad \delta \zstar
               = \delta z / r \quad \text{with} \quad r
-              = \frac{H + \eta}{H}
+              = \frac{H + \eta}{H} .
+\]
+Simple re-organisation of the above expressions gives
+\[
+  % \label{eq:PE_zstar_2}
+  \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) . 
 \]
 Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar,
@@ -776,15 +791,7 @@
 Also the divergence of the flow field is no longer zero as shown by the continuity equation:
 \[
-  \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) (r \; w *) = 0
+  \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 .
 \]
-
-% from MOM4p1 documentation
-To overcome problems with vanishing surface and/or bottom cells, we consider the zstar coordinate 
-\[
-  % \label{eq:PE_}
-  \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt)
-\]
-
-This coordinate is closely related to the "eta" coordinate used in many atmospheric models
+This \zstar coordinate is closely related to the "eta" coordinate used in many atmospheric models
 (see Black (1994) for a review of eta coordinate atmospheric models).
 It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves,
@@ -798,5 +805,5 @@
 These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to
 terrain following sigma models discussed in \autoref{subsec:PE_sco}.
-Additionally, since \zstar when $\eta = 0$,
+Additionally, since $\zstar = z$ when $\eta = 0$,
 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography.
 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of
@@ -804,21 +811,21 @@
 depending on the sophistication of the pressure gradient solver.
 The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of
-neutral physics parameterizations in \zstar models using the same techniques as in $z$-models
+neutral physics parameterizations in \zstar  models using the same techniques as in $z$-models
 (see Chapters 13-16 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$-models,
 as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO).
 
-The range over which \zstar varies is time independent $-H \leq \zstar \leq 0$.
-Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > ?H$.
+The range over which \zstar  varies is time independent $-H \leq \zstar \leq 0$.
+Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > -H$.
 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$.
 
-Because \zstar has a time independent range, all grid cells have static increments ds,
-and the sum of the ver tical increments yields the time independent ocean depth. %k ds = H.
+Because \zstar  has a time independent range, all grid cells have static increments ds,
+and the sum of the vertical increments yields the time independent ocean depth. %k ds = H.
 The \zstar coordinate is therefore invisible to undulations of the free surface,
 since it moves along with the free surface.
-This proper ty means that no spurious vertical transport is induced across surfaces of constant \zstar by
+This property means that no spurious vertical transport is induced across surfaces of constant \zstar  by
 the motion of external gravity waves.
-Such spurious transpor t can be a problem in z-models, especially those with tidal forcing.
-Quite generally, the time independent range for the \zstar coordinate is a very convenient property that
-allows for a nearly arbitrary ver tical resolution even in the presence of large amplitude fluctuations of
+Such spurious transport can be a problem in z-models, especially those with tidal forcing.
+Quite generally, the time independent range for the \zstar  coordinate is a very convenient property that
+allows for a nearly arbitrary vertical resolution even in the presence of large amplitude fluctuations of
 the surface height, again so long as $\eta > -H$.
 %end MOM doc %%%
@@ -870,9 +877,9 @@
 \begin{equation}
   \label{eq:PE_p_sco}
-  \nabla p |_z = \nabla p |_s - \pd[p]{s} \nabla z |_s
+  \nabla p |_z = \nabla p |_s - \frac{1}{e_3} \pd[p]{s} \nabla z |_s
 \end{equation}
 
 The second term in \autoref{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.
+leads to 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{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of the magnitude of this truncation error.
@@ -887,5 +894,5 @@
 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
+an 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 \autoref{sec:DOM_zgr}.
@@ -922,5 +929,5 @@
 
 The \ztilde -coordinate has been developed by \citet{leclair.madec_OM11}.
-It is available in \NEMO since the version 3.4.
+It is available in \NEMO since the version 3.4 and is more robust in version 4.0 than previously. 
 Nevertheless, it is currently not robust enough to be used in all possible configurations.
 Its use is therefore not recommended.
@@ -934,5 +941,5 @@
 \label{sec:PE_zdf_ldf}
 
-The primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than
+The hydrostatic primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than
 a few kilometres in the horizontal, a few meters in the vertical and a few minutes.
 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model.
@@ -984,5 +991,5 @@
 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
-(\eg Richardson number, Brunt-Vais\"{a}l\"{a} frequency...),
+(\eg Richardson number, Brunt-Vais\"{a}l\"{a} frequency, distance from the boundary ...),
 or computed from a turbulent closure model.
 The choices available in \NEMO are discussed in \autoref{chap:ZDF}).
@@ -1016,5 +1023,5 @@
 both horizontal and isoneutral operators have no effect on mean (\ie large scale) potential energy whereas
 potential energy is a main source of turbulence (through baroclinic instabilities).
-\citet{gent.mcwilliams_JPO90} have proposed a parameterisation of mesoscale eddy-induced turbulence which
+\citet{gent.mcwilliams_JPO90} proposed a 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.
@@ -1033,10 +1040,10 @@
 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.
+one uses an advective scheme which is diffusive enough to maintain the model stability.
 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of
 the advection scheme.
 
 All these parameterisations of subgrid scale physics have advantages and drawbacks.
-There are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are:
+They are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are:
 Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces,
 \citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes.
@@ -1141,5 +1148,5 @@
                          - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\
                       &= \lt(   \frac{1}{e_1}     \pd[ \lt( A^{lm}    \chi      \rt) ]{i} \rt.
-                              - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j}
+                              - \frac{1}{e_2 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{j} , 
                                 \frac{1}{e_2}     \pd[ \lt( A^{lm}    \chi      \rt) ]{j}
                          \lt. + \frac{1}{e_1 e_3} \pd[ \lt( A^{lm} \; e_3 \zeta \rt) ]{i} \rt)
@@ -1164,5 +1171,5 @@
 a geographical coordinate system \citep{lengaigne.madec.ea_JGR03}.
 
-\subsubsection{lateral bilaplacian momentum diffusive operator}
+\subsubsection{Lateral bilaplacian momentum diffusive operator}
 
 As for tracers, the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with
