Changeset 1224 for trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

Timestamp:

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

Author:

gm

Message:

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

File:

• trunk/DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (41 diffs)

trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

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