Index: branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2255)
+++ branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2268)
@@ -48,11 +48,10 @@
 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.
 
@@ -88,13 +87,13 @@
 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
+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.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -113,9 +112,9 @@
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
-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: 
@@ -167,4 +166,5 @@
 that are pure numerical artefacts. 
 
+\sgacomment{not sure whether 3 is still relevant after TRA-TRC branch is merged in?}
 % -------------------------------------------------------------------------------------------------------------
 %        2nd order centred scheme  
@@ -188,13 +188,14 @@
 (\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}
 
 % -------------------------------------------------------------------------------------------------------------
@@ -206,5 +207,5 @@
 
 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}
@@ -257,5 +258,5 @@
 \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{Zalesak_JCP79}. $c_u$ only departs from $1$ when the advective term 
@@ -264,8 +265,8 @@
 This scheme is tested and compared with MUSCL and the MPDATA scheme in 
 \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 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}$ 
+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 
@@ -302,5 +303,5 @@
 (\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.
 
@@ -339,6 +340,6 @@
 \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) 
+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 
@@ -352,4 +353,13 @@
 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 :
 
@@ -361,13 +371,5 @@
 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{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.
-
-(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\{	 
@@ -391,7 +393,6 @@
 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}
@@ -412,9 +413,9 @@
 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.
 
 
@@ -427,8 +428,9 @@
 
 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.
 
 % ================================================================
@@ -469,18 +471,18 @@
 \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.
+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. 
+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. 
+\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 
@@ -541,6 +543,5 @@
 background horizontal diffusion \citep{Guilyardi_al_CD01}. An alternative scheme 
 developed by \cite{Griffies_al_JPO98} which preserves both tracer and its variance 
-is currently been tested in \NEMO. It should be available in a forthcoming
-release.
+is also available in \NEMO (\np{ln\_traldf_grif}=true).
 
 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), the horizontal 
@@ -572,5 +573,5 @@
 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.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -584,13 +585,13 @@
 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. 
+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.
 
@@ -662,36 +663,38 @@
 the thickness of the top model layer. 
 
-Due to interactions and mass exchange with other media ($i.e.$ atmosphere, sea-ice, land),
+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 and not linked with $F_{mass}$, the water 
-exchange with the other media, and to the heat and salt content of this water exchange.
-In a forcoming release, these two parts, computed in the surface module (SBC), will included directly
+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.
-\gmcomment{  The specification of those fluxes is further detail in the SBC chapter (see \S\ref{SBC}).  }
-This change will provide a same forcing formulation for any tracers (including temperature and salinity).
+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. 
+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 cross the sea surface 
-(difference between the total surface heat flux and the fraction of the short wave flux that 
+$\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 exchanged
+$\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 evaporation-precipitation+freezing-melting because 
-the sea-ice is not currently embedded in the ocean but levitates above it. There is not mass
+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 link to the fact that sea-ice has a non-sero salinity, and the concentration/dilution effect
+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 
-a equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess, 
+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, \jp{lk\_vvl}=true):
+In the nonlinear free surface case (\key{vvl} is defined):
 \begin{equation} \label{Eq_tra_sbc}
 \begin{aligned}
@@ -704,5 +707,5 @@
 \end{equation} 
 
-In the linear free surface case (\key{vvl} not defined, , \jp{lk\_vvl}=false):
+In the linear free surface case (\key{vvl} not defined):
 \begin{equation} \label{Eq_tra_sbc_lin}
 \begin{aligned}
@@ -713,20 +716,19 @@
  \end{aligned}
 \end{equation} 
-where $\overline{x }^t$ means that $x$ is averaged over two consecutive time step 
-($t-\rdt/2$ and $t+\rdt/2$). Such a time averaged prevents the excitation of the 
+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 while their salinity is zero. Therefore, the heat content
-of \textit{emp} budget must be added to the temperature equation in variable volume case, 
-while it does not appear in constant volume. Similarly, the \textit{emp} budget affects 
-the ocean surface salinity in constant volume case (through the concentration dilution effect)
-while it does not appears explicitly in variable volume as salinity change will be
-induced by volume change. In both constant and variable volume, surface salinity 
-will change with ice-ocean salt flux and F/M flux without mass exchanges 
-($\textit{emp}_S - \textit{emp}$).
-
-Note that concentration/dilution effect due to F/M is computed using
+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} 
@@ -734,17 +736,17 @@
 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, since sea-surface salinity and F/M flux are 
+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.
+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 unbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$.
+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 that $(\partial_t\eta - \textit{emp})$ 
-and \textit{SSS} are correlated \citep{Roullet_Madec_JGR00}. 
-Note that, while quite small, the unbalance in constant volume case is larger 
-than the unbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}. 
-This is the reason why the modified filter is not applied in constant volume case.
+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.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -759,8 +761,8 @@
 
 When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true), 
-the solar radiation penetrates the top few 10 meters of the ocean, otherwise 
-all the heat flux is absorbed in the first ocean level (\np{ln\_flxqsr}=false). 
+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:
@@ -780,21 +782,21 @@
 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 10 centimetres. The fraction of $Q_{sr}$ 
+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 
-following a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
-of a few 10 centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the namtra\_qsr namlist).
+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 depths where it contributes to a penetrating flux of solar energy and thus 
-to local heating below the surface. 
-The way this second part of the solar energy penetrates in the ocean depends on 
-which formulation is chosen. In the simple 2-wavebands light penetration  (\np{ln\_qsr\_2bd}=true) 
-a chlorophyll-independent monochromatic formulation is also chosen for the shorter wavelengths, 
+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_0} + \left( 1-R\right) e^{-z / \xi_1} \right]
 \end{equation}
-where $\xi_1$ is the second extinction length scales associated with the shorter wavebands.  
-It is usually chosen to be 23~m through \np{rn\_si0} namelist parameter. 
+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).
@@ -802,26 +804,26 @@
 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 the 
-particules concentration and it is spectrally selective. \cite{Morel_JGR88} has shown 
+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 splitted into three wavebands: blue (400-500 nm), 
-green (500-600 nm) and red (600-700nm). For each wave-band, the chlorophyll-dependant 
+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 expression. As shown on Fig.\ref{Fig_traqsr_irradiance}, 
-this formulation, called RGB (Reed-Green-Blue), reproduces quite closely 
-the light penetration profiles predicted by the full spectal model with much faster 
-computing efficiently, in contrast with the 2-bands formulation. 
+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 inverse of the extinction length scales) are tabulated over 61 nonuniform 
+($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 type of chlorophyll can be used in the RGB formulation:
-(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) observed 
-time varying chlorophyll (\np{nn\_chdta}=0) ; (3) simulated time varying chlorophyll
-by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the later case, the RGB 
-formulation is used to calculated both the phytoplankton light limitation in PISCES 
+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. 
 
@@ -829,8 +831,8 @@
 is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
 
-When $z$-coordinate is preferred to $s$-coordinate, the depth of $w-$levels does 
+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 until that level. 
+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 
@@ -841,7 +843,7 @@
 \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 wavebands chlorophyll-independant formulation (blue), 
-a chlorophyll-dependant monochromatic formulation (green), 4 waveband RGB formulation (red), 
+\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}.}
@@ -873,10 +875,10 @@
 earth cooling. This flux is weak compared to surface fluxes (a mean 
 global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), but it is 
-systematically positive and acts on the densest water masses. 
+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_OS09}. 
 
-The presence or not of geothermal heating is controlled by the namelist 
+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 
@@ -904,9 +906,9 @@
 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 
+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 in large 
-proportions downstream of a sill \citep{Willebrand_al_PO01}, and the thickness 
+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. 
 
@@ -914,9 +916,9 @@
 \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 
-one of the improvment introduced by \citet{Campin_Goosse_Tel99}.
+the improvements introduced by \citet{Campin_Goosse_Tel99}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -927,5 +929,5 @@
 
 When applying sigma-diffusion (\key{trabbl} defined and \np{nn\_bbl\_ldf} set to 1), 
-the diffusive flux between two adjacent cells living at the ocean bottom is given by 
+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
@@ -933,5 +935,5 @@
 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 
 and  $A_l^\sigma$ the lateral diffusivity in the BBL. Following \citet{Beckmann_Doscher1997}, 
-the latter is prescribed with a spatial dependence, $e.g.$ in the conditional form
+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}
@@ -943,6 +945,6 @@
 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist 
 parameter \np{rn\_ahtbbl} and usually set to a value much larger 
-than the one used on lateral mixing in open ocean. The constraint in \eqref{Eq_tra_bbl_coef} 
-implies that sigma-like diffusion only occurs when density above the sea floor, at the top of 
+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, 
@@ -951,6 +953,6 @@
 	\nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S
 \end{equation} 
-where $\rho$, $\alpha$ and $\beta$ are function of $\overline{T}^\sigma$, 
-$\overline{S}^\sigma$, $\overline{H}^\sigma$, the along bottom mean temperature, 
+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.
 
@@ -961,14 +963,16 @@
 \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.7\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
-\caption{Advective/diffusive Bottom Boundary Layer. The BBl parameterisation is 
+\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 downsloping flow is defined either as the transport of the bottom 
+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 connecting directly $kup$ and $kdwn$
+The green arrow indicates the diffusive BBL flux directly connecting $kup$ and $kdwn$
 ocean bottom cells.
 connection}
@@ -988,5 +992,5 @@
 water to move down the slope. 
 
-\np{nn\_bbl\_adv} = 1 : the downsloping velocity is chosen to be the Eulerian
+\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
@@ -995,30 +999,30 @@
 greater depth ($i.e.$ $\vect{U}  \cdot  \nabla H>0$).
 
-\np{nn\_bbl\_adv} = 2 : the downsloping velocity is chosen to be proportional to $\Delta \rho$,
-the density difference between the top and down cell densities \citep{Campin_Goosse_Tel99}.
+\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 
-downsloping flow, here in the $i$-direction (Fig.\ref{Fig_bbl}), is simply given by 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 second, is the coefficient of proportionality 
+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 top and bottom cells, respectively.
+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 retained. The possible values for $\gamma$ 
+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}.  
 
-The scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ 
+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 downsloping plume and the surrounding 
-water at intermediate depth. The entrainment is replaced by the vertical mixing 
-included in the advection scheme. Let us consider as an example the 
-case display in Fig.\ref{Fig_bbl} where the density at level $(i,kup)$ is 
+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 downsloping flow \eqref{Eq_bbl_dw}, 
+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: 
@@ -1062,7 +1066,7 @@
 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}. 
 
@@ -1074,10 +1078,10 @@
 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 
@@ -1085,5 +1089,5 @@
 \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 $S_o$ are used 
+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 
@@ -1091,8 +1095,8 @@
 structure in equilibrium with its physics. 
 The choice of the shape of the Newtonian damping is controlled by two 
-namelist parameters \np{nn\_zdmp}. The former allows to specified the 
-width of the equatorial band in which no damping is applied as well as a decrease 
-in vicinity of the coast and a damping everywhere in the Red and Med Seas, 
-whereas the latter set a damping acting in the mixed layer or not. 
+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 
