Changeset 2282 for branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex

Timestamp:

2010-10-15T16:42:00+02:00 (14 years ago)

Author:

gm

Message:

ticket:#658 merge DOC of all the branches that form the v3.3 beta

File:

• branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (58 diffs)

branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_TRA.tex

@@ -8 +8 @@
 % missing/update 
 % traqsr: need to coordinate with SBC module
-% trabbl : advective case to be discussed
-%        diffusive case : add : only the bottom ocean cell is concerned
-%        ==> addfigure on bbl
 
 %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
@@ -51 +48 @@
 In the present chapter we also describe the diagnostic equations used to compute 
 the sea-water properties (density, Brunt-Vais\"{a}l\"{a} frequency, specific heat and 
-freezing point) although the associated modules ($i.e.$ \mdl{eosbn2}, \mdl{ocfzpt} 
-and \mdl{phycst}) are (temporarily) located in the NEMO/OPA directory.
-
-The different options available to the user are managed by namelist logical or 
+freezing point with associated modules \mdl{eosbn2}, \mdl{ocfzpt} and \mdl{phycst}).
+
+The different options available to the user are managed by namelist logicals or 
 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 
+where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 
+The CPP key (when it exists) is \textbf{key\_trattt}. The equivalent code can be 
 found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
 
@@ -70 +66 @@
       {Tracer Advection (\mdl{traadv})}
 \label{TRA_adv}
-%------------------------------------------nam_traadv-----------------------------------------------------
-\namdisplay{nam_traadv}
+%------------------------------------------namtra_adv-----------------------------------------------------
+\namdisplay{namtra_adv}
 %-------------------------------------------------------------------------------------------------------------
 
@@ -77 +73 @@
 fluxes. Its discrete expression is given by :
 \begin{equation} \label{Eq_tra_adv}
-ADV_\tau =-\frac{1}{b_T} \left( 
+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. 
+-\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 
-\left[ {e_{2u} {\kern 1pt}{\kern 1pt}\;u\;\tau _u } \right]+\delta _j \left[ 
-{e_{1v} {\kern 1pt}v\;\tau _v } \right]\;} \right)-\frac{1}{\mathop 
-e\nolimits_{3T} }\delta _k \left[ {w\;\tau _w } \right]
+ADV_\tau = - \frac{1}{e_{1t}\,e_{2t}} \left( \; \delta_i \left[ e_{2u} \;u \;\tau_u \right] 
+                                                                 + \delta_j \left[ e_{1v} \;v \;\tau_v  \right] \; \right)
+                   -  \frac{1}{e_{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. Indeed, it is obtained
+$e_{3u} =e_{3v} =e_{3t} $. The flux form in \eqref{Eq_tra_adv} 
+implicitly requires 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. 
+$\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant volume (default option) 
+or variable volume (\key{vvl} defined) 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 
+by replacing $\tau$ by the number 1 in (\ref{Eq_tra_adv}) we recover the discrete form of 
 the continuity equation which is used to calculate the vertical velocity.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -117 +112 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-The key difference between the advection schemes used in \NEMO is the choice 
+The key difference between the advection schemes available in \NEMO is the choice 
 made in space and time interpolation to define the value of the tracer at the 
 velocity points (Fig.~\ref{Fig_adv_scheme}). 
 
-Along solid lateral and bottom boundaries a zero tracer flux is naturally 
+Along solid lateral and bottom boundaries a zero tracer flux is automatically 
 specified, since the normal velocity is zero there. At the sea surface the 
 boundary condition depends on the type of sea surface chosen: 
 \begin{description}
-\item  [rigid-lid formulation:] $w=0$ at the surface, so the advective 
-fluxes through the surface are zero.
 \item [linear free surface:] the first level thickness is constant in time: 
 the vertical boundary condition is applied at the fixed surface $z=0$ 
 rather than on the moving surface $z=\eta$. There is a non-zero advective 
-flux which is set for all advection schemes as the product of surface 
-velocity (at $z=0$) by the first level tracer value: 
-$\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $. 
+flux which is set for all advection schemes as 
+$\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$ 
+the product of surface velocity (at $z=0$) by the first level tracer value.
 \item [non-linear free surface:] (\key{vvl} is defined) 
 convergence/divergence in the first ocean level moves the free surface 
@@ -144 +137 @@
 term is the product of the time derivative of the tracer and the free surface 
 height, two quantities that are not correlated (see \S\ref{PE_free_surface}, 
-and also \citet{Roullet2000,Griffies2001,Campin2004}).
+and also \citet{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}).
 
 The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco}) 
 is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}). 
-When advective bottom boundary layer (\textit{bbl}) and/or eddy induced velocity 
-(\textit{eiv}) parameterisations are used it is the \textit{now} \textit{effective} 
-velocity ($i.e.$ the sum of the eulerian, the bbl and/or the eiv velocities) which is used.
-
-The choice of an advection scheme is made in the \np{nam\_traadv} namelist, by 
+When eddy induced velocity (\textit{eiv}) parameterisation is used it is the \textit{now} 
+\textit{effective} velocity ($i.e.$ the sum of the eulerian and eiv velocities) which is used.
+
+The choice of an advection scheme is made in the \textit{nam\_traadv} namelist, by 
 setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The 
 corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 
@@ -174 +166 @@
 that are pure numerical artefacts. 
 
+\sgacomment{not sure whether 3 is still relevant after TRA-TRC branch is merged in?}
 % -------------------------------------------------------------------------------------------------------------
 %        2nd order centred scheme  
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})]
-         {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=.true.)}
+         {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=true)}
 \label{TRA_adv_cen2}
 
@@ -195 +188 @@
 (\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 
+it is advantageous to combine the \textit{cen2} scheme with an upstream scheme
+in specific areas which require 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).
+This combined scheme has been included for specific grid points in the ORCA2 and ORCA4 configurations only.
 
 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. Note also that 
+have this order of accuracy. \gmcomment{Note also that ... blah, blah}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -209 +203 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})]
-           {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=.true.)}
+           {$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=true)}
 \label{TRA_adv_cen4}
 
 In the $4^{th}$ order formulation (to be implemented), tracer values are 
-evaluated at velocity points as a $4^{th}$ order interpolation, and thus uses 
+evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on 
 the four neighbouring $T$-points. For example, in the $i$-direction:
 \begin{equation} \label{Eq_tra_adv_cen4}
@@ -247 +241 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]
-         {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=.true.)}
+         {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}
 \label{TRA_adv_tvd}
 
@@ -264 +258 @@
 \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 
+There exist many ways to define $c_u$, each corresponding 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 
+\citet{Zalesak_JCP79}. $c_u$ only departs from $1$ when the advective term 
 produces a local extremum in the tracer field. The resulting scheme is quite 
 expensive but \emph{positive}. It can be used on both active and passive tracers. 
 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. 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}) 
-$\tau _u^{cen2}$ is evaluated using the \textit{now} tracer while $\tau _u^{ups}$ 
+\citet{Levy_al_GRL01}; note that in this paper it is referred to as "FCT" (Flux corrected 
+transport) rather than TVD. The TVD scheme is implemented in the \mdl{traadv\_tvd} module.
+
+For stability reasons (see \S\ref{DOM_nxt}),
+$\tau _u^{cen2}$ is evaluated  in (\ref{Eq_tra_adv_tvd}) using the \textit{now} tracer while $\tau _u^{ups}$ 
 is evaluated using the \textit{before} tracer. In other words, the advective part of 
 the scheme is time stepped with a leap-frog scheme while a forward scheme is 
@@ -287 +281 @@
 
 The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been 
-implemented by \citet{Levy2001}. In its formulation, the tracer at velocity points 
+implemented by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points 
 is evaluated assuming a linear tracer variation between two $T$-points 
 (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_muscl}
    \tau _u^{mus} = \left\{      \begin{aligned}
-         &\tau _i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\Delta t}{e_{1u}} \right)
+         &\tau _i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
          &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\
-         &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\Delta t}{e_{1u} } \right)
+         &\tau _{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right)
          &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0
    \end{aligned}    \right.
@@ -306 +300 @@
 For an ocean grid point adjacent to land and where the ocean velocity is 
 directed toward land, two choices are available: an upstream flux 
-(\np{ln\_traadv\_muscl}=.true.) or a second order flux 
-(\np{ln\_traadv\_muscl2}=.true.). Note that the latter choice does not ensure 
+(\np{ln\_traadv\_muscl}=true) or a second order flux 
+(\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. The two MUSCL schemes are computed 
+on both active and passive tracers. The two MUSCL schemes are implemented 
 in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules.
 
@@ -316 +310 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs})]
-         {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=.true.)}
+         {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=true)}
 \label{TRA_adv_ubs}
 
@@ -333 +327 @@
 
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation 
-error \citep{Sacha2005}. The overall performance of the advection 
+error \citep{Shchepetkin_McWilliams_OM05}. 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. 
@@ -344 +338 @@
 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..
-
-For stability reasons  (see \S\ref{DOM_nxt}), in \eqref{Eq_tra_adv_ubs}, 
-the first term (which corresponds to a second order centred scheme) 
+\np{ln\_traadv\_ubs}=true.
+
+For stability reasons  (see \S\ref{DOM_nxt}),
+the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order centred scheme) 
 is evaluated using the \textit{now} tracer (centred in time) while the 
 second term (which is the diffusive part of the scheme), is 
 evaluated using the \textit{before} tracer (forward in time). 
-This choice is discussed by \citet{Webb1998} in the context of the 
+This choice is discussed by \citet{Webb_al_JAOT98} 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}. 
+leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
 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.
+
+Four different options are possible for the vertical 
+component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated 
+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{Shchepetkin_McWilliams_OM05} 
+implementation of UBS in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case 
+has dispersion properties similar to an eighth-order accurate conventional scheme.
+The current reference version uses method b)
 
 Note that :
@@ -368 +371 @@
 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)} 
-a TVD scheme, or  \textit{(c)} an interpolation based on conservative 
-parabolic splines following the \citet{Sacha2005} implementation of UBS 
-in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case has dispersion properties 
-similar to an eighth-order accurate conventional scheme.
-
-(3) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
+(2) 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\{  
@@ -398 +393 @@
 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy 
 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_traadv_ubs2}. This should be 
- changed in forthcoming release.
+ $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that NEMO v3.3 still uses 
+ \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}.
  %%%
  \gmcomment{the change in UBS scheme has to be done}
@@ -409 +403 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [QUICKEST scheme (QCK) (\np{ln\_traadv\_qck})]
-         {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=.true.)}
+         {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=true)}
 \label{TRA_adv_qck}
 
@@ -419 +413 @@
 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 
+However, 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.
+This no longer guarantees the positivity of the scheme. The use of TVD in the vertical 
+direction (as for the UBS case) should be implemented to restore this property.
 
 
@@ -430 +424 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})]
-         {Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=.true.)}
+         {Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=true)}
 \label{TRA_adv_ppm}
 
 The Piecewise Parabolic Method (PPM) proposed by Colella and Woodward (1984) 
-is based on a quadradic piecewise rebuilding. Like the QCK scheme, it is associated 
+\sgacomment{reference?}
+is based on a quadradic piecewise construction. Like the QCK scheme, it is associated 
 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 3.0.
+version 3.3.
 
 % ================================================================
@@ -446 +441 @@
 \label{TRA_ldf}
 %-----------------------------------------nam_traldf------------------------------------------------------
-\namdisplay{nam_traldf}
+\namdisplay{namtra_ldf}
 %-------------------------------------------------------------------------------------------------------------
  
@@ -465 +460 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap})]
-         {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=.true.) }
+         {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=true) }
 \label{TRA_ldf_lap}
 
@@ -471 +466 @@
 surfaces is given by: 
 \begin{equation} \label{Eq_tra_ldf_lap}
-D_T^{lT} =\frac{1}{b_T} \left( \;
+D_T^{lT} =\frac{1}{b_tT} \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.
+where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells. 
+It is implemented 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 \np{ln\_traldf\_hor}=\np{ln\_zco}=.true.. 
+or without partial steps, 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 \np{ln\_traldf\_hor}=\np{ln\_zco}=true. 
 In both cases, it significantly contributes to diapycnal mixing. 
 It is therefore not recommended.
 
 Note that 
-(a) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}=e_{3v}=e_{3T}$, 
+(a) In the 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 
+(b) In the 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 
@@ -499 +494 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})]
-         {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=.true.)}
+         {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=true)}
 \label{TRA_ldf_iso}
 
@@ -507 +502 @@
 \begin{equation} \label{Eq_tra_ldf_iso}
 \begin{split}
- 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]
+ 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.    \\ 
@@ -525 +520 @@
  \end{split}
  \end{equation}
-where $b_T= e_{1T}\,e_{2T}\,e_{3T}$  is the volume of $T$-cells, 
+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., 
-or both \np{ln\_traldf\_hor}=.true. and \np{ln\_zco}=.true.. The way these 
+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 
 and lateral boundaries, the turbulent fluxes of heat and salt are set to zero 
@@ -546 +541 @@
 of the tracer variance. Nevertheless the treatment performed on the slopes 
 (see \S\ref{LDF}) allows the model to run safely without any additional 
-background horizontal diffusion \citep{Guily2001}. An alternative scheme 
-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 
+background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme 
+developed by \cite{Griffies_al_JPO98} which preserves both tracer and its variance 
+is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of 
+the algorithm is given in App.\ref{Apdx_Griffies}.
+
+Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal 
 derivatives at the bottom level in \eqref{Eq_tra_ldf_iso} require a specific 
 treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
@@ -559 +554 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap})]
-         {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=.true.)}
+         {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=true)}
 \label{TRA_ldf_bilap}
 
 The lateral fourth order bilaplacian operator on tracers is obtained by 
-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 
-\np{ln\_zco}=.false.. In both cases, it can contribute diapycnal mixing, 
+applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption 
+on boundary conditions: both 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 
+\np{ln\_zco}=false. In both cases, it can contribute diapycnal mixing, 
 although less than in the laplacian case. It is therefore not recommended.
 
@@ -579 +574 @@
 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.
+number of code-lines.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -585 +580 @@
 % -------------------------------------------------------------------------------------------------------------
 \subsection   [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})]
-         {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=.true.)}
+         {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=true)}
 \label{TRA_ldf_bilapg}
 
@@ -591 +586 @@
 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 can be found in the
+coast, normal to the bottom and normal to 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. 
+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. 
+This rotated bilaplacian operator has never been seriously 
+tested. There are no guarantees that it is either free of bugs or correctly formulated. 
 Moreover, the stability range of such an operator will be probably quite 
-narrow, requiring a significantly smaller time-step than the one used on 
+narrow, requiring a significantly smaller time-step than the one used with an
 unrotated operator.
 
@@ -618 +613 @@
 \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}
@@ -641 +636 @@
 The large eddy coefficient found in the mixed layer together with high 
 vertical resolution implies that in the case of explicit time stepping 
-(\np{ln\_zdfexp}=.true.) there would be too restrictive a constraint on 
+(\np{ln\_zdfexp}=true) there would be too restrictive a constraint on 
 the time step. Therefore, the default implicit time stepping is preferred 
 for the vertical diffusion since it overcomes the stability constraint. 
-A forward time differencing scheme (\np{ln\_zdfexp}=.true.) using a time 
-splitting technique (\np{n\_zdfexp} $> 1$) is provided as an alternative. 
-Namelist variables \np{ln\_zdfexp} and \np{n\_zdfexp} apply to both 
+A forward time differencing scheme (\np{ln\_zdfexp}=true) using a time 
+splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 
+Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both 
 tracers and dynamics. 
 
@@ -667 +662 @@
 enhance readability of the code. The two formulations are completely 
 equivalent; the forcing terms in trasbc are the surface fluxes divided by 
-the thickness of the top model layer. Following \citet{Roullet2000} the 
-forcing on an ocean tracer, $c$, can be split into two parts: $F_{ext}$, the 
-flux of tracer crossing the sea surface and not linked with the water 
-exchange with the atmosphere, $F_{wf}^d$, and $F_{wf}^i$ the forcing 
-on the tracer concentration associated with this water exchange. 
-The latter forcing has two components: a direct effect of change 
-in concentration associated with the tracer carried by the water flux, 
-and an indirect concentration/dilution effect :
-\begin{equation*}
-\begin{split}
- F^C &= F_{ext} + F_{wf}^d                                          +F_{wf}^i    \\
-        &= F_{ext} - \left( c_E \, E - c_p \,P - c_R \,R \right) +c\left( E-P-R \right)
-\end{split}
-\end{equation*} 
-
-\gmcomment{add here a description of the variable names used in the above equation}
-
-Two cases must be distinguished, the nonlinear free surface case 
-(\key{vvl} is defined) and the linear free surface case. The first case 
-is simpler, because the indirect concentration/dilution effect is naturally 
-taken into account by letting the vertical scale factors vary in time. 
-The salinity of water exchanged at the surface is assumed to be zero, 
-so there is no salt flux at the free surface, except in the presence of sea ice. 
-The heat flux at the free surface is the sum of $F_{ext}$, the direct 
-heating/cooling (by the total non-penetrative heat flux) and $F_{wf}^e$ 
-the heat carried by the water exchanged through the surface (evaporation, 
-precipitation, runoff). The temperature of precipitation is not well known. 
-In the model we assume that this water has the same temperature as 
-the sea surface temperature. The resulting forcing terms for temperature 
-T and salinity S are: 
-\begin{equation} \label{Eq_tra_forcing}
+the thickness of the top model layer. 
+
+Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, sea-ice, land),
+the change in the heat and salt content of the surface layer of the ocean is due both 
+to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
+ and to the heat and salt content of the mass exchange.
+\sgacomment{ the following does not apply to the release to which this documentation is 
+attached and so should not be included ....
+In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly
+in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux.
+The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}). 
+This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity).
+ 
+In the current version, the situation is a little bit more complicated. }
+The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following 
+forcing fields (used on tracers):
+
+$\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 
+(i.e. the difference between the total surface heat flux and the fraction of the short wave flux that 
+penetrates into the water column, see \S\ref{TRA_qsr})
+
+$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
+
+$\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchange
+
+$\bullet$ \textit{rnf}, the mass flux associated with runoff (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+
+The $\textit{emp}_S$ field is not simply the budget of evaporation-precipitation+freezing-melting because 
+the sea-ice is not currently embedded in the ocean but levitates above it. There is no mass
+exchanged between the sea-ice and the ocean. Instead we only take into account the salt
+flux associated with the non-zero salinity of sea-ice, and the concentration/dilution effect
+due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into 
+an equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess, 
+the surface boundary condition on temperature and salinity is applied as follows:
+
+In the nonlinear free surface case (\key{vvl} is defined):
+\begin{equation} \label{Eq_tra_sbc}
 \begin{aligned}
- &F^T =\frac{Q_{ns} }{\rho _o \;C_p \,e_{3T} }-\frac{\text{EMP}\;\left. T 
-\right|_{k=1} }{e_{3T} }  & \\ 
-\\
-& F^S =\frac{\text{EMP}_S\;\left. S \right|_{k=1} }{e_{3T} }   &
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }   
+           &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t  & \\ 
+%
+& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 
+           &\overline{ \left( (\textit{emp}_S - \textit{emp})\;\left. S \right|_{k=1}  \right) }^t   & \\   
  \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 
-non-penetrative part of the net surface heat flux (difference between 
-the total surface heat flux and the fraction of the short wave flux that 
-penetrates into the water column), the product EMP$_S\;.\left. S \right|_{k=1}$ 
-is  the ice-ocean salt flux, and $\left. S\right|_{k=1}$ is the sea surface 
-salinity (\textit{SSS}). The total salt content is conserved in this formulation 
-(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?     ===> gm :  NO to be added at NOCS 
-
-In the second case (linear free surface), the vertical scale factors are 
-fixed in time so that the concentration/dilution effect must be added in 
-the \mdl{trasbc} module. Because of the hypothesis made for the 
-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} }
+
+In the linear free surface case (\key{vvl} not defined):
+\begin{equation} \label{Eq_tra_sbc_lin}
+\begin{aligned}
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns} }^t  & \\ 
+%
+& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} } 
+           &\overline{ \left( \textit{emp}_S\;\left. S \right|_{k=1}  \right) }^t   & \\   
+ \end{aligned}
 \end{equation} 
-
-The salinity forcing is still given by \eqref{Eq_tra_forcing} but the 
-definition of EMP$_S$ is different: it is the total surface freshwater 
-budget (evaporation minus precipitation minus river runoff plus 
-the rate of change of the sea ice thickness). The total salt content 
-is not exactly conserved (\citet{Roullet2000}. 
-See also \S\ref{PE_free_surface}).
-
-In the case of the rigid lid approximation, the surface salinity forcing $F^s$ 
-is also expressed by \eqref{Eq_tra_forcing}, but now the global integral of 
-the product of EMP and S, is not compensated by the advection of fluid 
-through the top level: this is because in the rigid lid case \textit{w(k=1) = 0} 
-(in contrast to the linear free surface case). As a result, even if the budget 
-of \textit{EMP} is zero on average over the whole ocean domain, the 
-associated salt flux is not, since sea-surface salinity and \textit{EMP} are 
-intrinsically correlated (high \textit{SSS} are found where evaporation is 
-strong whilst low \textit{SSS} is usually associated with high precipitation 
-or river runoff).
-
-The $Q_{ns} $ and \textit{EMP} fields are defined and updated in the 
-\mdl{sbcmod} module (see \S\ref{SBC}).
+where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps 
+($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the 
+divergence of odd and even time step (see \S\ref{STP}).
+
+The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained 
+by assuming that the temperature of precipitation and evaporation are equal to
+the ocean surface temperature and that their salinity is zero. Therefore, the heat content
+of the \textit{emp} budget must be added to the temperature equation in the variable volume case, 
+while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects 
+the ocean surface salinity in the constant volume case (through the concentration dilution effect)
+while it does not appears explicitly in the variable volume case since salinity change will be
+induced by volume change. In both constant and variable volume cases, surface salinity 
+will change with ice-ocean salt flux and F/M flux (both contained in $\textit{emp}_S - \textit{emp}$) without mass exchanges.
+
+Note that the concentration/dilution effect due to F/M is computed using
+a constant ice salinity as well as a constant ocean salinity. 
+This approximation suppresses the correlation between \textit{SSS} 
+and F/M flux, allowing the ice-ocean salt exchanges to be conservative.
+Indeed, if this approximation is not made, even if the F/M budget is zero 
+on average over the whole ocean domain and over the seasonal cycle, 
+the associated salt flux is not zero, since sea-surface salinity and F/M flux are 
+intrinsically correlated (high \textit{SSS} are found where freezing is 
+strong whilst low \textit{SSS} is usually associated with high melting areas).
+
+Even using this approximation, an exact conservation of heat and salt content 
+is only achieved in the variable volume case. In the constant volume case, 
+there is a small imbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$.
+Nevertheless, the salt content variation is quite small and will not induce
+a long term drift as there is no physical reason for $(\partial_t\eta - \textit{emp})$ 
+and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}. 
+Note that, while quite small, the imbalance in the constant volume case is larger 
+than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 
+This is the reason why the modified filter is not applied in the constant volume case.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -753 +758 @@
 \label{TRA_qsr}
 %--------------------------------------------namqsr--------------------------------------------------------
-\namdisplay{namqsr}
+\namdisplay{namtra_qsr}
 %--------------------------------------------------------------------------------------------------------------
 
-When the penetrative solar radiation option is used (\np{ln\_flxqsr}=.true.), 
-the solar radiation penetrates the top few meters of the ocean, otherwise 
-all the heat flux is absorbed in the first ocean level (\np{ln\_flxqsr}=.false.). 
+When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true), 
+the solar radiation penetrates the top few tens of meters of the ocean. If it is not used 
+(\np{ln\_flxqsr}=false) all the heat flux is absorbed in the first ocean level. 
 Thus, in the former case a term is added to the time evolution equation of 
-temperature \eqref{Eq_PE_tra_T} whilst the surface boundary condition is 
+temperature \eqref{Eq_PE_tra_T} and the surface boundary condition is 
 modified to take into account only the non-penetrative part of the surface 
 heat flux:
@@ -769 +774 @@
 \end{split}
 \end{equation}
-
-where $I$ is the downward irradiance. The additional term in \eqref{Eq_PE_qsr} 
-is discretized as follows:
+where $Q_{sr}$ is the penetrative part of the surface heat flux ($i.e.$ the shortwave radiation) 
+and $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$). 
+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]
-\end{equation}
-
-A formulation involving two extinction coefficients is assumed for the 
-downward irradiance $I$ \citep{Paulson1977}:
+\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]
+\end{equation}
+
+The shortwave radiation,  $Q_{sr}$, consists of energy distributed across a wide spectral range. 
+The ocean is strongly absorbing for wavelengths longer than 700~nm and these 
+wavelengths contribute to heating the upper few tens of centimetres. The fraction of $Q_{sr}$ 
+that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ (specified 
+through namelist parameter \np{rn\_abs}).  It is assumed to penetrate the ocean 
+with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
+of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namelist).
+For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy 
+propagates to larger depths where it contributes to 
+local heating. 
+The way this second part of the solar energy penetrates into the ocean depends on 
+which formulation is chosen. In the simple 2-waveband light penetration scheme  (\np{ln\_qsr\_2bd}=true) 
+a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 
+leading to the following expression  \citep{Paulson1977}:
 \begin{equation} \label{Eq_traqsr_iradiance}
-I(z) = Q_{sr} \left[Re^{-z / \xi_1} + \left( 1-R\right) e^{-z / \xi_2} \right]
-\end{equation}
-where $Q_{sr}$ is the penetrative part of the surface heat flux, 
-$\xi_1$ and $\xi_2$ are two extinction length scales and $R$ 
-determines the relative contribution of the two terms. 
-The default values used correspond to a Type I water in Jerlov's [1968] 
-%
-\gmcomment : Jerlov reference to be added
-%
-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), 
-so all the solar radiation that reaches the last ocean level is absorbed 
-in that level. The trend in \eqref{Eq_tra_qsr} associated with the 
-penetration of the solar radiation is added to the temperature trend, 
-and the surface heat flux is modified in routine \mdl{traqsr}. 
-Note that in the $z$-coordinate, the depth of $T-$levels depends 
-on the single variable $k$. A one dimensional array of the coefficients 
-$gdsr(k) = Re^{-z_w (k)/\xi_1} + (1-R)e^{-z_w (k)/\xi_2}$ can then 
-be computed once and saved in memory. Moreover \textit{nksr}, 
-the level at which $gdrs$ becomes negligible (less than the 
-computer precision) is computed once, and the trend associated 
-with the penetration of the solar radiation is only added until that level. 
-Finally, note that when the ocean is shallow (< 200~m), part of the 
+I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
+\end{equation}
+where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths.  
+It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 
+The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in 
+Jerlov's (1968) classification (oligotrophic waters).
+
+Such assumptions have been shown to provide a very crude and simplistic 
+representation of observed light penetration profiles (\cite{Morel_JGR88}, see also 
+Fig.\ref{Fig_traqsr_irradiance}). Light absorption in the ocean depends on 
+particle concentration and is spectrally selective. \cite{Morel_JGR88} has shown 
+that an accurate representation of light penetration can be provided by a 61 waveband 
+formulation. Unfortunately, such a model is very computationally expensive. 
+Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this 
+formulation in which visible light is split into three wavebands: blue (400-500 nm), 
+green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependent 
+attenuation coefficient is fitted to the coefficients computed from the full spectral model 
+of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}), assuming 
+the same power-law relationship. As shown in Fig.\ref{Fig_traqsr_irradiance}, 
+this formulation, called RGB (Red-Green-Blue), reproduces quite closely 
+the light penetration profiles predicted by the full spectal model, but with much greater 
+computational efficiency. The 2-bands formulation does not reproduce the full model very well. 
+
+The RGB formulation is used when \np{ln\_qsr\_rgb}=true. The RGB attenuation coefficients
+($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform 
+chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
+in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation:
+(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed 
+time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll
+by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB 
+formulation is used to calculate both the phytoplankton light limitation in PISCES 
+or LOBSTER and the oceanic heating rate. 
+
+The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation 
+is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
+
+When the $z$-coordinate is preferred to the $s$-coordinate, the depth of $w-$levels does 
+not significantly vary with location. The level at which the light has been totally 
+absorbed ($i.e.$ it is less than the computer precision) is computed once, 
+and the trend associated with the penetration of the solar radiation is only added down to that level. 
+Finally, note that when the ocean is shallow ($<$ 200~m), part of the 
 solar radiation can reach the ocean floor. In this case, we have 
 chosen that all remaining radiation is absorbed in the last ocean 
-level ($i.e.$ $I_w$ is masked). 
-
-When coupling with a biological model (for example PISCES or LOBSTER), 
-it is possible to calculate the light attenuation using information from 
-the biology model. Without biological model, it is still possible to introduce 
-a horizontal variation of the light attenuation by using the observed ocean 
-surface color. At the time of writing, the latter has not been implemented
- in the reference version.
-%
-\gmcomment{  {yellow}{case 4 bands and bio-coupling to add !!!}  }
-%
+level ($i.e.$ $I$ is masked). 
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] \label{Fig_traqsr_irradiance}  \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_Irradiance.pdf}
+\caption{Penetration profile of the downward solar irradiance 
+calculated by four models. Two waveband chlorophyll-independent formulation (blue), 
+a chlorophyll-dependent monochromatic formulation (green), 4 waveband RGB formulation (red), 
+61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 
+(a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. From \citet{Lengaigne_al_CD07}.}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % -------------------------------------------------------------------------------------------------------------
@@ -824 +859 @@
 \label{TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
-\namdisplay{nambbc}
+\namdisplay{namtra_bbc}
 %--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \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}$) used by \cite{Emile-Geay_Madec_OSD08}.
-It is 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_OS09}.
+It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -840 +875 @@
 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. 
+global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it is 
+systematically positive \sgacomment{positive definite?} 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 
-parameter  \np{ngeo\_flux}. If this parameter is set to 1, a constant 
+Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}. 
+
+The presence of geothermal heating is controlled by the namelist 
+parameter  \np{nn\_geoflx}. When this parameter is set to 1, a constant 
 geothermal heating is introduced whose value is given by the 
-\np{ngeo\_flux\_const}, which is also a namelist parameter. If it is set to 2, 
+\np{nn\_geoflx\_cst}, which is also a namelist parameter. When it is set to 2, 
 a spatially varying geothermal heat flux is introduced which is provided 
-in the geothermal\_heating.nc NetCDF file (Fig.\ref{Fig_geothermal}).
+in the \ifile{geothermal\_heating} NetCDF file (Fig.\ref{Fig_geothermal}).
 
 % ================================================================
 % Bottom Boundary Layer
 % ================================================================
-\section  [Bottom Boundary Layer (\textit{trabbl}, \textit{trabbl\_adv} )]
-      {Bottom Boundary Layer (\mdl{trabbl}, \mdl{trabbl\_adv})}
+\section  [Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})]
+      {Bottom Boundary Layer (\mdl{trabbl} - \key{trabbl})}
 \label{TRA_bbl}
 %--------------------------------------------nambbl---------------------------------------------------------
@@ -865 +900 @@
 In a $z$-coordinate configuration, the bottom topography is represented by a 
 series of discrete steps. This is not adequate to represent gravity driven 
-downslope flows. Such flows arise downstream of sills such as the Strait of 
-Gibraltar, Bab El Mandeb, or Denmark Strait, where dense water formed in 
-marginal seas flows into a basin filled with less dense water. The amount of 
-entrainment that occurs in these gravity plumes is critical in determining the 
-density and volume flux of the densest waters of the ocean, such as 
-Antarctic Bottom Water, or North Atlantic Deep Water. $z$-coordinate 
-models tend to overestimate the entrainment, because the gravity flow is 
-mixed down vertically by convection as it goes ``downstairs'' following the 
-step topography, sometimes over a thickness much larger than the thickness 
-of the observed gravity plume. A similar problem occurs in the $s$-coordinate when 
-the thickness of the bottom level varies in large proportions downstream of 
-a sill \citep{Willebrand2001}, and the thickness of the plume is not resolved. 
-
-The idea of the bottom boundary layer (BBL) parameterisation first introduced by 
-\citet{BeckDos1998} is to allow a direct communication between 
+downslope flows. Such flows arise either downstream of sills such as the Strait of 
+Gibraltar or Denmark Strait, where dense water formed in marginal seas flows 
+into a basin filled with less dense water, or along the continental slope when dense 
+water masses are formed on a continental shelf. The amount of entrainment 
+that occurs in these gravity plumes is critical in determining the density 
+and volume flux of the densest waters of the ocean, such as Antarctic Bottom Water, 
+or North Atlantic Deep Water. $z$-coordinate models tend to overestimate the 
+entrainment, because the gravity flow is mixed vertically by convection 
+as it goes ''downstairs'' following the step topography, sometimes over a thickness 
+much larger than the thickness of the observed gravity plume. A similar problem 
+occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly 
+downstream of a sill \citep{Willebrand_al_PO01}, and the thickness 
+of the plume is not resolved. 
+
+The idea of the bottom boundary layer (BBL) parameterisation, first introduced by 
+\citet{Beckmann_Doscher1997}, is to allow a direct communication between 
 two adjacent bottom cells at different levels, whenever the densest water is 
-located above the less dense water. The communication can be by a diffusive 
-(diffusive BBL), advective fluxes (advective BBL), or both. In the current 
+located above the less dense water. The communication can be by a diffusive flux
+(diffusive BBL), an advective flux (advective BBL), or both. In the current 
 implementation of the BBL, only the tracers are modified, not the velocities. 
 Furthermore, it only connects ocean bottom cells, and therefore does not include 
-the improvment proposed by \citet{Campin_Goosse_Tel99}.
+the improvements introduced by \citet{Campin_Goosse_Tel99}.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Diffusive Bottom Boundary layer (\key{bbl\_diff})}
+\subsection{Diffusive Bottom Boundary layer (\np{nn\_bbl\_ldf}=1)}
 \label{TRA_bbl_diff}
 
-When applying sigma-diffusion (\key{trabbl} is defined), the diffusive flux between 
-two adjacent cells living at the ocean bottom is given by 
+When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 
+the diffusive flux between two adjacent cells at the ocean floor is given by 
 \begin{equation} \label{Eq_tra_bbl_diff}
 {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T
 \end{equation} 
 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 
-and  $A_l^\sigma $ the lateral diffusivity in the BBL. Following \citet{BeckDos1998}, 
-the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form
+and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997}, 
+the latter is prescribed with a spatial dependence, $i.e.$ in the conditional form
 \begin{equation} \label{Eq_tra_bbl_coef}
 A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
@@ -909 +945 @@
 \end{equation} 
 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist 
-parameter \np{atrbbl}. $A_{bbl}$ is usually set to a value much larger 
-than the one used on lateral mixing in open ocean. 
-Note that in practice, \eqref{Eq_tra_bbl_coef} constraint is applied 
-separately in the two horizontal directions, and the density gradient in 
-\eqref{Eq_tra_bbl_coef} is evaluated at $\overline{H}^i$ ($\overline{H}^j$) 
-using the along bottom mean temperature and salinity. 
+parameter \np{rn\_ahtbbl} and usually set to a value much larger 
+than the one used for lateral mixing in the open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
+implies that sigma-like diffusion only occurs when the density above the sea floor, at the top of 
+the slope, is larger than in the deeper ocean (see green arrow in Fig.\ref{Fig_bbl}). 
+In practice, this constraint is applied separately in the two horizontal directions, 
+and the density gradient in \eqref{Eq_tra_bbl_coef} is evaluated with the log gradient formulation: 
+\begin{equation} \label{Eq_tra_bbl_Drho}
+   \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
+\end{equation} 
+where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$, 
+$\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, 
+salinity and depth, respectively.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection   {Advective Bottom Boundary Layer (\key{bbl\_adv})}
+\subsection   {Advective Bottom Boundary Layer  (\np{nn\_bbl\_adv}= 1 or 2)}
 \label{TRA_bbl_adv}
 
+\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following
+if this is not what is meant then "downwards sloping flow" is also a possibility"}
 
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_bbl}  \begin{center}
-\includegraphics[width=0.8\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
-\caption{Advective Bottom Boundary Layer.}
+\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
+\caption{Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is 
+activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 
+Red arrows indicate the additional overturning circulation due to the advective BBL. 
+The transport of the downslope flow is defined either as the transport of the bottom 
+ocean cell (black arrow), or as a function of the along slope density gradient. 
+The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$
+ocean bottom cells.
+connection}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
+
+%!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
+%!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
+%!!        i.e. transport proportional to the along-slope density gradient
+
 %%%gmcomment   :  this section has to be really written
 
-The advective BBL is in fact not only an advective one but include a diffusive 
-component as we chose an upstream scheme to perform the advection within 
-the BBL. The associated diffusion only act in the stream direction and is 
-proportional to the velocity.
-
-When applying sigma-advection (\key{trabbl\_adv} defined), the advective 
-flux between two adjacent cells living at the ocean bottom is given by 
-\begin{equation} \label{Eq_tra_bbl_fadv}
-{\rm {\bf F}}_\sigma={\rm {\bf U}}_h^\sigma \; \overline{T}^\sigma
-\end{equation} 
-with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 
-and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{BeckDos1998}, 
-the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form
-\begin{equation} \label{Eq_tra_bbl_Aadv}
-A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
- A_{bbl} \quad \quad \mbox{if}     \quad    \nabla_\sigma \rho \cdot \nabla H<0 
-              \quad \quad \mbox{and} \quad         {\rm {\bf U}}_h  \cdot \nabla H<0 \\ 
- \\
- 0\quad \quad \;\,\mbox{otherwise} \\ 
- \end{array}} \right.
-\end{equation} 
+When applying an advective BBL (\np{nn\_bbl\_adv} = 1 or 2), an overturning 
+circulation is added which connects two adjacent bottom grid-points only if dense 
+water overlies less dense water on the slope. The density difference causes dense 
+water to move down the slope. 
+
+\np{nn\_bbl\_adv} = 1 : the downslope velocity is chosen to be the Eulerian
+ocean velocity just above the topographic step (see black arrow in Fig.\ref{Fig_bbl}) 
+\citep{Beckmann_Doscher1997}. It is a \textit{conditional advection}, that is, advection
+is allowed only if dense water overlies less dense water on the slope ($i.e.$ 
+$\nabla_\sigma \rho  \cdot  \nabla H<0$) and if the velocity is directed towards 
+greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$).
+
+\np{nn\_bbl\_adv} = 2 : the downslope velocity is chosen to be proportional to $\Delta \rho$,
+the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
+The advection is allowed only  if dense water overlies less dense water on the slope ($i.e.$ 
+$\nabla_\sigma \rho  \cdot  \nabla H<0$). For example, the resulting transport of the 
+downslope flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by the 
+following expression:
+\begin{equation} \label{Eq_bbl_Utr}
+ u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right)
+\end{equation}
+where $\gamma$, expressed in seconds, is the coefficient of proportionality 
+provided as \np{rn\_gambbl}, a namelist parameter, and \textit{kup} and \textit{kdwn} 
+are the vertical index of the higher and lower cells, respectively.
+The parameter $\gamma$ should take a different value for each bathymetric 
+step, but for simplicity, and because no direct estimation of this parameter is 
+available, a uniform value has been assumed. The possible values for $\gamma$ 
+range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}.  
+
+Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ 
+using the upwind scheme. Such a diffusive advective scheme has been chosen 
+to mimic the entrainment between the downslope plume and the surrounding 
+water at intermediate depths. The entrainment is replaced by the vertical mixing 
+implicit in the advection scheme. Let us consider as an example the 
+case displayed in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is 
+larger than the one at level $(i,kdwn)$. The advective BBL scheme
+modifies the tracer time tendency of the ocean cells near the 
+topographic step by the downslope flow \eqref{Eq_bbl_dw}, 
+the horizontal \eqref{Eq_bbl_hor}  and the upward \eqref{Eq_bbl_up} 
+return flows as follows: 
+\begin{align} 
+\partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw}
+                                     +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right)  \label{Eq_bbl_dw} \\
+%
+\partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
+               + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{Eq_bbl_hor} \\
+%
+\intertext{and for $k =kdw-1,\;..., \; kup$ :} 
+%
+\partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k}
+               + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{Eq_bbl_up}
+\end{align}
+where $b_t$ is the $T$-cell volume. 
+
+Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in 
+the model outputs. It has to be used to compute the effective velocity 
+as well as the effective overturning circulation.
 
 % ================================================================
@@ -960 +1052 @@
       {Tracer damping (\mdl{tradmp})}
 \label{TRA_dmp}
-%--------------------------------------------namtdp-----------------------------------------------------
-\namdisplay{namtdp}
+%--------------------------------------------namtra_dmp-------------------------------------------------
+\namdisplay{namtra_dmp}
 %--------------------------------------------------------------------------------------------------------------
 
@@ -969 +1061 @@
 \begin{split}
  \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right)  \\
-\\ 
  \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right) 
  \end{split}
@@ -976 +1067 @@
 are given temperature and salinity fields (usually a climatology). 
 The restoring term is added when \key{tradmp} is defined. 
-It also requires that both \key{temdta} and \key{saldta} are defined 
+It also requires that both \key{dtatem} and \key{dtasal} are defined 
 ($i.e.$ that $T_o$ and $S_o$ are read). The restoring coefficient 
-$S_o$ is a three-dimensional array initialized by the user in routine 
+$\gamma$ is a three-dimensional array initialized by the user in routine 
 \rou{dtacof} also located in module \mdl{tradmp}. 
 
@@ -988 +1079 @@
 field for a passive tracer study). The first case applies to regional 
 models that have artificial walls instead of open boundaries. 
-In the vicinity of these walls, $S_o$ takes large values (equivalent to 
+In the vicinity of these walls, $\gamma$ takes large values (equivalent to 
 a time scale of a few days) whereas it is zero in the interior of the 
 model domain. The second case corresponds to the use of the robust 
 diagnostic method \citep{Sarmiento1982}. It allows us to find the velocity 
-field consistent with the model dynamics whilst having a $T$-$S$ field 
-close to a given climatological field ($T_o -S_o$). The time scale 
+field consistent with the model dynamics whilst having a $T$, $S$ field 
+close to a given climatological field ($T_o$, $S_o$). The time scale 
 associated with $S_o$ is generally not a constant but spatially varying 
 in order to respect other properties. For example, it is usually set to zero 
 in the mixed layer (defined either on a density or $S_o$ criterion) 
-\citep{Madec1996} and in the equatorial region 
-\citep{Reverdin1991, Fujio1991, MartiTh1992} since these two regions 
-have a short time scale of adjustment; while smaller $S_o$ are used 
+\citep{Madec_al_JPO96} and in the equatorial region 
+\citep{Reverdin1991, Fujio1991, Marti_PhD92} since these two regions 
+have a short time scale of adjustment; while smaller $\gamma$ are used 
 in the deep ocean where the typical time scale is long \citep{Sarmiento1982}. 
 In addition the time scale is reduced (even to zero) along the western 
 boundary to allow the model to reconstruct its own western boundary 
-structure in equilibrium with its physics. The choice of a 
-Newtonian damping acting in the mixed layer or not is controlled by
-namelist parameter \np{nmldmp}.
+structure in equilibrium with its physics. 
+The choice of the shape of the Newtonian damping is controlled by two 
+namelist parameters \np{??} and \np{nn\_zdmp}. The former allows us to specify: the 
+width of the equatorial band in which no damping is applied; a decrease 
+in the vicinity of the coast; and a damping everywhere in the Red and Med Seas.
+The latter sets whether damping should act in the mixed layer or not. 
+The time scale associated with the damping depends on the depth as
+a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as 
+bottom value and a transition depth of \np{rn\_dep}.  
 
 The robust diagnostic method is very efficient in preventing temperature 
@@ -1013 +1110 @@
 by stabilising the water column too much.
 
-An example of the computation of $S_o$ for robust diagnostic experiments 
+An example of the computation of $\gamma$ for a robust diagnostic experiment 
 with the ORCA2 model is provided in the \mdl{tradmp} module 
 (subroutines \rou{dtacof} and \rou{cofdis} which compute the coefficient 
@@ -1029 +1126 @@
 %--------------------------------------------------------------------------------------------------------------
 
-The general framework for tracer time stepping is a leap-frog scheme, 
-$i.e.$ a three level centred time scheme associated with a Asselin time 
-filter (cf. \S\ref{DOM_nxt}):
+The general framework for tracer time stepping is a modified leap-frog scheme 
+\citep{Leclair_Madec_OM09}, $i.e.$ a three level centred time scheme associated 
+with a Asselin time filter (cf. \S\ref{STP_mLF}):
 \begin{equation} \label{Eq_tra_nxt}
-\begin{split}
-T^{t+\Delta t} &= T^{t-\Delta t} + 2 \, \Delta t  \ \text{RHS}_T^t   \\
+\begin{aligned}
+(e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t & \\
 \\
-T_f^t  \;\ \quad &= T^t \;\quad +\gamma \,\left[ {T_f^{t-\Delta t} -2T^t+T^{t+\Delta t}} \right]
-\end{split}
+(e_{3t}T)_f^t  \;\ \quad &= (e_{3t}T)^t \;\quad 
+                                    &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] &  \\
+                                 & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]  &                      
+\end{aligned}
 \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 
-coefficient. $\gamma$ is initialized as \np{atfp} (\textbf{namelist} parameter). 
-Its default value is \np{atfp=0.1}. 
+where RHS is the right hand side of the temperature equation, 
+the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient,
+and $S$ is the total forcing applied on $T$ ($i.e.$ fluxes plus content in mass exchanges). 
+$\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 
+Its default value is \np{rn\_atfp}=$10^{-3}$. Note that the forcing correction term in the filter
+is not applied in linear free surface (\jp{lk\_vvl}=false) (see \S\ref{TRA_sbc}.
+Not also that in constant volume case, the time stepping is performed on $T$, 
+not on its content, $e_{3t}T$.
 
 When the vertical mixing is solved implicitly, the update of the \textit{next} tracer 
@@ -1049 +1152 @@
 
 In order to prepare for the computation of the \textit{next} time step, 
-a swap of tracer arrays is performed: $T^{t-\Delta t} = T^t$ and $T^t = T_f$. 
+a swap of tracer arrays is performed: $T^{t-\rdt} = T^t$ and $T^t = T_f$. 
 
 % ================================================================
@@ -1064 +1167 @@
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Equation of State (\np{neos} = 0, 1 or 2)}
+\subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)}
 \label{TRA_eos}
 
 It is necessary to know the equation of state for the ocean very accurately 
 to determine stability properties (especially the Brunt-Vais\"{a}l\"{a} frequency), 
-particularly in the deep ocean. The ocean density is a non linear empirical 
-function of \textit{in situ }temperature, salinity and pressure. The reference 
-equation of state is that defined by the Joint Panel on Oceanographic Tables 
-and Standards \citep{UNESCO1983}. It was the standard equation of state 
-used in early releases of OPA. However, even though this computation is 
-fully vectorised, it is quite time consuming ($15$ to $20${\%} of the total 
-CPU time) since it requires the prior computation of the \textit{in situ} 
-temperature from the model \textit{potential} temperature using the 
-\citep{Bryden1973} polynomial for adiabatic lapse rate and a $4^th$ order 
-Runge-Kutta integration scheme. Since OPA6, we have used the 
-\citet{JackMcD1995} equation of state for seawater instead. It allows the 
-computation of the \textit{in situ} ocean density directly as a function of 
-\textit{potential} temperature relative to the surface (an \NEMO variable), 
-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 
-the values of the various coefficients have been adjusted by \citet{JackMcD1995} 
-in order to directly use the \textit{potential} temperature instead of the 
-\textit{in situ} one. This reduces the CPU time of the in situ density computation 
-to about $3${\%} of the total CPU time, while maintaining a quite accurate 
-equation of state.
-
-In the computer code, a \textit{true} density $d$ is computed, $i.e.$ the ratio 
-of seawater volumic mass to $\rho_o$, a reference volumic mass (\textit{rau0} 
-defined in \mdl{phycst}, usually $rau0= 1,020~Kg/m^3$). The default option 
-(namelist prameter \np{neos}=0) is the \citet{JackMcD1995} equation of state. 
-Its use is highly recommended. However, for process studies, it is often 
-convenient to use a linear approximation of the density$^{\ast}$
-\footnote{$^{\ast }$ With the linear equation of state there is no longer 
-a distinction between \textit{in situ} and \textit{potential} density. Cabling 
-and thermobaric effects are also removed.}. 
-Two linear formulations are available: a function of $T$ only (\np{neos}=1) 
-and a function of both $T$ and $S$ (\np{neos}=2):
+particularly in the deep ocean. The ocean seawater volumic mass, $\rho$, 
+abusively called density, is a non linear empirical function of \textit{in situ} 
+temperature, salinity and pressure. The reference equation of state is that 
+defined by the Joint Panel on Oceanographic Tables and Standards 
+\citep{UNESCO1983}. It was the standard equation of state used in early 
+releases of OPA. However, even though this computation is fully vectorised, 
+it is quite time consuming ($15$ to $20${\%} of the total CPU time) since 
+it requires the prior computation of the \textit{in situ} temperature from the 
+model \textit{potential} temperature using the \citep{Bryden1973} polynomial 
+for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme. 
+Since OPA6, we have used the \citet{JackMcD1995} equation of state for 
+seawater instead. It allows the computation of the \textit{in situ} ocean density 
+directly as a function of \textit{potential} temperature relative to the surface 
+(an \NEMO variable), 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 the values of the various coefficients have 
+been adjusted by \citet{JackMcD1995} in order to directly use the \textit{potential} 
+temperature instead of the \textit{in situ} one. This reduces the CPU time of the 
+\textit{in situ} density computation to about $3${\%} of the total CPU time, 
+while maintaining a quite accurate equation of state.
+
+In the computer code, a \textit{true} density anomaly, $d_a= \rho / \rho_o - 1$, 
+is computed, with $\rho_o$ a reference volumic mass. Called \textit{rau0} 
+in the code, $\rho_o$ is defined in \mdl{phycst}, and a value of $1,035~Kg/m^3$. 
+This is a sensible choice for the reference density used in a Boussinesq ocean 
+climate model, as, with the exception of only a small percentage of the ocean, 
+density in the World Ocean varies by no more than 2$\%$ from $1,035~kg/m^3$ 
+\citep{Gill1982}.
+
+The default option (namelist parameter \np{nn\_eos}=0) is the \citet{JackMcD1995} 
+equation of state. Its use is highly recommended. However, for process studies, 
+it is often convenient to use a linear approximation of the density.
+With such an equation of state there is no longer a distinction between 
+\textit{in situ} and \textit{potential} density and both cabbeling and thermobaric
+effects are removed.
+Two linear formulations are available: a function of $T$ only (\np{nn\_eos}=1) 
+and a function of both $T$ and $S$ (\np{nn\_eos}=2):
 \begin{equation} \label{Eq_tra_eos_linear}
-\begin{aligned}
- d(T)    &= {\rho (T)} / {\rho _0 } &&= 1.028 - \alpha \;T     \\ 
- d(T,S) &= {\rho (T,S)}                &&= \ \ \ \beta \;S - \alpha \;T 
-\end{aligned}
+\begin{split}
+  d_a(T)       &=  \rho (T)      /  \rho_o   - 1     =  \  0.0285         -  \alpha   \;T     \\ 
+  d_a(T,S)    &=  \rho (T,S)   /  \rho_o   - 1     =  \  \beta \; S       -  \alpha   \;T    
+\end{split}
 \end{equation} 
 where $\alpha$ and $\beta$ are the thermal and haline expansion 
 coefficients, and $\rho_o$, the reference volumic mass, $rau0$. 
-($\alpha$ and $\beta$ can be modified through the \np{ralpha} and 
-\np{rbeta} namelist parameters). Note that when $d$ is a function 
-of $T$ only (\np{neos}=1), the salinity is a passive tracer and can be 
+($\alpha$ and $\beta$ can be modified through the \np{rn\_alpha} and 
+\np{rn\_beta} namelist parameters). Note that when $d_a$ is a function 
+of $T$ only (\np{nn\_eos}=1), the salinity is a passive tracer and can be 
 used as such.
 
@@ -1119 +1227 @@
 %        Brunt-Vais\"{a}l\"{a} Frequency
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{neos} = 0, 1 or 2)}
+\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
 \label{TRA_bn2}
 
@@ -1128 +1236 @@
  iso-neutral diffusion). In particular, one must be aware that $N^2$ has to 
  be computed with an \textit{in situ} reference. The expression for $N^2$ 
- depends on the type of equation of state used (\np{neos} namelist parameter).
-
-For \np{neos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} 
+ depends on the type of equation of state used (\np{nn\_eos} namelist parameter).
+
+For \np{nn\_eos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987} 
 polynomial expression is used (with the pressure in decibar approximated by 
 the depth in meters): 
@@ -1137 +1245 @@
       \left(  \alpha / \beta \ \delta_{k+1/2}[T]     - \delta_{k+1/2}[S]   \right) 
 \end{equation} 
-where $\alpha$ ($\beta$) is the thermal (haline) expansion coefficient. 
-They are a function of  
-$\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$, 
-and  $z_w$, with $T$ the \textit{potential} temperature and 
-$\widetilde{S}$ a salinity anomaly. 
+where $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
+They are a function of  $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$, 
+and  $z_w$, with $T$ the \textit{potential} temperature and $\widetilde{S}$ a salinity anomaly. 
 Note that both $\alpha$ and $\beta$ depend on \textit{potential} 
 temperature and salinity which are averaged at $w$-points prior 
@@ -1147 +1253 @@
 then averaged to $w$-points.
 
-When a linear equation of state is used (\np{neos}=1 or 2, 
+When a linear equation of state is used (\np{nn\_eos}=1 or 2, 
 \eqref{Eq_tra_bn2} reduces to:
 \begin{equation} \label{Eq_tra_bn2_linear}
@@ -1158 +1264 @@
 %        Specific Heat
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Specific Heat (\textit{phycst})]
+\subsection    [Specific Heat (\textit{phycst})]
          {Specific Heat (\mdl{phycst})}
 \label{TRA_adv_ldf}
@@ -1171 +1277 @@
 Its value is set in \mdl{phycst} module. 
 
-%%%
-\gmcomment{ STEVEN:  consistency, no other computer variable names are 
-supplied, so why this one}
-%%%
 
 % -------------------------------------------------------------------------------------------------------------
 %        Freezing Point of Seawater
 % -------------------------------------------------------------------------------------------------------------
-\subsection   [Freezing Point of Seawater (\textit{ocfzpt})]
-         {Freezing Point of Seawater (\mdl{ocfzpt})}
+\subsection   [Freezing Point of Seawater]
+         {Freezing Point of Seawater}
 \label{TRA_fzp}
 
@@ -1194 +1296 @@
 \eqref{Eq_tra_eos_fzp} is only used to compute the potential freezing point of 
 sea water ($i.e.$ referenced to the surface $p=0$), thus the pressure dependent 
-terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The \textit{before} 
-and \textit{now} surface freezing point is introduced in the code as $fzptb$ and 
-$fzptn$ 2D arrays together with a  \textit{now} mask (\textit{freezn}) which takes 
-the value 0 or 1 depending on whether the ocean temperature is above or at the 
-freezing point. Caution: do not confuse \textit{freezn} with the fraction of lead 
-(\textit{frld}) defined in LIM.  
-
-%%%
-\gmcomment{STEVEN: consistency, not many computer variable names are supplied, so why these    ===>  gm  I agree   this should evolve both here and in the code itself}
-%%%
+terms in \eqref{Eq_tra_eos_fzp} (last term) have been dropped. The freezing
+point is computed through \textit{tfreez}, a \textsc{Fortran} function that can be found 
+in \mdl{eosbn2}.  
 
 % ================================================================
@@ -1214 +1309 @@
 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial bottom cells (\np{ln\_zps}=.true.), in general, tracers in horizontally 
+With partial bottom cells (\np{ln\_zps}=true), in general, tracers in horizontally 
 adjacent cells live at different depths. Horizontal gradients of tracers are needed 
 for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure 
@@ -1228 +1323 @@
 \begin{figure}[!p] \label{Fig_Partial_step_scheme}  \begin{center}
 \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Partial_step_scheme.pdf}
-\caption{ Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate (\np{ln\_zps}=.true.) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$.  }
+\caption{ Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate (\np{ln\_zps}=true) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. The horizontal difference is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$.  }
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>