Index: /branches/DEV_r1826_DOC/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- /branches/DEV_r1826_DOC/DOC/TexFiles/Biblio/Biblio.bib	(revision 2194)
+++ /branches/DEV_r1826_DOC/DOC/TexFiles/Biblio/Biblio.bib	(revision 2195)
@@ -412,5 +412,5 @@
 }
 
-@ARTICLE{Bernie_al_OM05,
+@ARTICLE{Bernie_al_JC05,
   author = {D. J. Bernie and S. J. Woolnough and J. M. Slingo and E. Guilyardi},
   title = {Modelling diurnal and intraseasonal variability of the ocean mixed
@@ -2624,4 +2624,25 @@
   owner = {gm},
   timestamp = {2008.02.09}
+}
+
+@ARTICLE{Morel_JGR88,
+  author = {A. Morel},
+  title = {Optical modeling of the upper ocean in relation to its biogenous
+	matter content (Case I waters)},
+  journal = JGR,
+  year = {1988},
+  volume = {93},
+  pages = {10,749--10,768},
+  owner = {gm},
+  timestamp = {2010.10.07}
+}
+
+@ARTICLE{Morel_Maritorena_JGR01,
+  author = {A. Morel and S. Maritorena},
+  title = {Bio-optical properties of oceanic waters: a reappraisal},
+  journal = JGR,
+  year = {2001},
+  owner = {gm},
+  timestamp = {2010.10.07}
 }
 
Index: /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2194)
+++ /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 2195)
@@ -1131,5 +1131,5 @@
 The general framework for dynamics time stepping is a leap-frog scheme, 
 $i.e.$ a three level centred time scheme associated with an Asselin time filter 
-(cf. \S\ref{DOM_nxt}). The scheme is applied to the velocity except when using 
+(cf. Chap.\ref{STP}). The scheme is applied to the velocity except when using 
 the flux form of momentum advection (cf. \S\ref{DYN_adv_cor_flux}) in variable 
 volume level case (\key{vvl} defined), where it has to be applied to the thickness 
@@ -1154,5 +1154,7 @@
 where RHS is the right hand side of the momentum equation, the subscript $f$ 
 denotes filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is 
-initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = 0.1.
+initialized as \np{nn\_atfp} (namelist parameter). Its default value is \np{nn\_atfp} = $10^{-3}$.
+In both cases, the modified Asselin filter is not applied since a perfect conservation 
+is not an issue for momentum equation.
 
 Note that with the filtered free surface, the update of the \textit{next} velocities 
Index: /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2194)
+++ /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 2195)
@@ -103,11 +103,6 @@
 In the real ocean, $\textit{emp}=\textit{emp}_S$ and the ocean salt content is conserved, 
 but it exist several numerical reasons why this equality should be broken. 
-For example:
-
-When the rigid-lid assumption is made, the ocean volume becomes constant and 
-thus, $\textit{emp}=0$, not $\textit{emp}_S$.
-
-When the ocean is coupled to a sea-ice model, the water exchanged between ice and 
-ocean is slightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case, 
+For example, when the ocean is coupled to a sea-ice model, the water exchanged between 
+ice and ocean is slightly salty (mean sea-ice salinity is $\sim $\textit{4 psu}). In this case, 
 $\textit{emp}_{S}$ take into account both concentration/dilution effect associated with 
 freezing/melting and the salt flux between ice and ocean, while \textit{emp} is 
@@ -359,5 +354,5 @@
 
 % ================================================================
-%        Addition of river runoffs
+%        River runoffs
 % ================================================================
 \section   [river runoffs (\textit{sbcrnf})]
@@ -375,15 +370,44 @@
 coastal modelling and becomes more and more often open ocean and climate modelling 
 \footnote{At least a top cells thickness of 1~meter and a 3 hours forcing frequency are
-required to properly represent the diurnal cycle \citep{Bernie_al_OM05}.}.
-
-
- 
+required to properly represent the diurnal cycle \citep{Bernie_al_JC05}. see also \S\ref{SBC_dcy}.}.
+
+
+To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the \mdl{tra\_sbc} module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below:
+
+
+Rachel:
+
 It is convenient to introduce the river runoff in the model as a surface 
-fresh water flux. 
-
-
-%Griffies:  River runoff generally enters the ocean at a nonzero depth rather than through the surface. Many global models, however, have traditionally inserted river runoff to the top model cell. Such can become problematic numerically and physically when the top grid cells are reÞned to levels common in coastal modelling. Hence, more applications are now considering the input of runoff throughout a nonzero depth. Likewise, sea ice can melt at depth, thus necessitating a mass transport to occur within the ocean between the liquid and solid water masses.
-
-\colorbox{yellow}{Nevertheless, Pb of vertical resolution and increase of Kz in vicinity of }
+fresh water flux. This is the defualt option within NEMO, and there is then
+ the option for the user to increase vertical mixing in the vicinity of the rivermouth.
+
+However, this method is not very appropriate for coastal modelling.  As such its now possible
+ to specify, in a netcdf input file, the temperature and salinity of the river, along with the 
+depth (in metres) which the river should be added to.  This enables to river to be correctly 
+added through the water column, instead of as a surface flux, and also means the temperature 
+or salinity (for low salinity outflow) of the river impacts the surrounding ocean.
+
+For temperature -999 is taken as missing data and the river temperature is taken to be the 
+surface temperatue at the river point.  For the depth parameter a value of -1 means the 
+river is added to the surface box only, and a value of -999 means the river is added through 
+the entire water column.
+
+Namelist options, \np{ln\_rnf\_depth}, \np{ln\_rnf\_sal} and \np{ln\_rnf\_temp} control whether 
+the river attributes (depth, salinity and temperature) are read in and used.  If these are set 
+as false the river is added to the surface box only, assumed to be fresh (0~psu), and/or 
+taken as surface temperature respectively.
+
+It is also possible for runnoff to be specified as a negative value for modelling flow through 
+straits, i.e. modelling the baltic flow in and out of the North Sea.  When the flow is out of the 
+domain there is no change in temperature and salinity, regardless of the namelist options used.
+
+The runoff value and attributes are read in in sbcrnf.  The mass/volume addition is added to the 
+divergence term in \rou{sbc\_rnf\_div}.  The dilution effect of the river is automatically applied through 
+the vertical tracer advection, and the direct flux of tracers into the domain is done in trasbc.
+
+
+\colorbox{yellow}{Nevertheless, Pb of vertical resolution and 3D input : increase vertical mixing near river mouths to mimic a 3D river 
+
+All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface.}
 
 \colorbox{yellow}{river mouths{\ldots}}
@@ -396,9 +420,51 @@
 
 In the current \NEMO setup river runoff is added to emp fluxes, these are then applied at just the sea surface as a volume change (in the variable volume case this is a literal volume change, and in the linear free surface case the free surface is moved) and a salt flux due to the concentration/dilution effect.  There is also an option to increase vertical mixing near river mouths; this gives the effect of having a 3d river.  All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface.
-Our aim was to code the option to specify the temperature and salinity of river runoff, (as well as the amount), along with the depth that the river water will affect.  This would make it possible to model low salinity outflow, such as the Baltic, and would allow the ocean temperature to be affected by river runoff.  The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between.
+Our aim was to code the option to specify the temperature and salinity of river runoff, (as well as the amount), along with the depth that the river water will affect.  This would make it possible to model low salinity outflow, such as the Baltic, and would allow the ocean temperature to be affected by river runoff.  
+
+The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between.
 
 To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the tra_sbc module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below:
 
 }
+
+
+% ================================================================
+%        Diurnal cycle
+% ================================================================
+\section   [Diurnal  cycle (\textit{sbcdcy})]
+			{Diurnal cycle (\mdl{sbcdcy})}
+\label{SBC_dcy}
+%------------------------------------------namsbc_rnf----------------------------------------------------
+%\namdisplay{namsbc} 
+%-------------------------------------------------------------------------------------------------------------
+
+\cite{Bernie_al_JC05} have shown that to capture 90$\%$ of the diurnal variability of 
+SST requires a vertical resolution in upper ocean of 1~m or better and a temporal resolution 
+of the surface fluxes of 3~h or less. Unfortunately high frequency forcing fields are rare, 
+not to say inexistent. Nevertheless, it is possible to obtain a reasonable diurnal cycle 
+of the SST knowning only short wave flux (SWF) at high frequency \citep{Bernie_al_CD07}.
+Furthermore, only the knowledge of daily mean value of SWF is needed, 
+as higher frequency variations can be reconstructed from them, assuming that 
+the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle 
+of incident SWF. The \cite{Bernie_al_CD07} reconstruction algorithm is available
+in \NEMO by setting \np{ln\_dm2dc}=true (a \textit{namsbc} namelist parameter) when using 
+CORE bulk formulea (\np{ln\_blk\_core}=true) or the flux formulation (\np{ln\_flx}=true). 
+The algorithm used is detailed in the appendix~A of \cite{Bernie_al_CD07} 
+and illustrated on Fig.\ref{Fig_SBC_diurnal}. 
+
+Note that the reconstruction can change the daily mean value of SWF by a few tenth of 
+W/m$^2$ ($\pm 0.5$~W/m$^2$ with a 1~h sampling) but this is not a systematic bias. 
+Note also that the setting a diurnal cycle in SWF is highly recommended  when 
+the top layer thickness approach 1~m or less, otherwise large error in SST can 
+appear due to an inconsistency between the scale of the vertical resolution 
+and the forcing acting on that scale.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] \label{Fig_SBC_diurnal}  \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_SBC_diurnal.pdf}
+\caption{Example of recontruction of the diurnal cycle variation of short wave flux  
+from daily mean values. From \citet{Bernie_al_CD07}.}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % ================================================================
Index: /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2194)
+++ /branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 2195)
@@ -666,61 +666,85 @@
 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.
-These two parts are included in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux.
-The specification of those fluxes is further detail in the SBC chapter (see \S\ref{SBC}).
-
-In the nonlinear free surface case (\key{vvl} is defined), the forcing term of $T$ and $S$
-conservation equations is simply given by:
+In a forcoming release, these two parts, computed in the surface module (SBC), will 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).
+ 
+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 
+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{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
+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
+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, 
+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):
 \begin{equation} \label{Eq_tra_sbc}
 \begin{aligned}
- F^T &=\frac{ Q_{ns} }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }      \\ 
+ &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{ F_{salt} }{\rho _o \; \left. e_{3t} \right|_{k=1} }   
+& 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 both $Q_{ns}$ and $F_{salt}$ are averaged over two consecutive time step (before and now)
-in order to prevent the exitation of the divergence of odd and even time step (see \S\ref{STP}).
-
-In the linear free surface case (\key{vvl} is not defined), the ocean volume is assumed 
-to be fixed in time, and so are the vertical scale factors. In order to compensate the 
-volume change due to water exchange with the other media, an extra term must be added,
-the so-called concentration/dilution effect. The forcing term becomes:
-\begin{equation} \label{Eq_tra_sbc_linssh}
+
+In the linear free surface case (\key{vvl} not defined, , \jp{lk\_vvl}=false):
+\begin{equation} \label{Eq_tra_sbc_lin}
 \begin{aligned}
- &F^T = \frac{ \overline{Q_{ns}}^t                                     }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }   
-          -  \frac{ \overline{F_{mass}}^t\;\left. T \right|_{k=1} }{\rho _o \; \left. e_{3t} \right|_{k=1} }   \\ 
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns} }^t  & \\ 
 %
-& F^S =\frac{ \overline{F_{salt}}^t                                     }{\rho _o \,\left. e_{3t} \right|_{k=1} } 
-          -  \frac{ \overline{F_{mass}}^t\;\left. S \right|_{k=1} }{\rho _o \; \left. e_{3t} \right|_{k=1} }  \\   
+& 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 total salt content is no more exactly conserved (\citet{Roullet_Madec_JGR00}. 
-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 \textit{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 \textit{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 $\textit{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 
+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 
+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
+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, 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 unbalance 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.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -735,5 +759,5 @@
 
 When the penetrative solar radiation option is used (\np{ln\_flxqsr}=true), 
-the solar radiation penetrates the top few meters of the ocean, otherwise 
+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). 
 Thus, in the former case a term is added to the time evolution equation of 
@@ -747,51 +771,81 @@
 \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}:
+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}$ 
+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).
+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, 
+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{rn\_si1}, \np{rn\_si2} and \np{rn\_abs} 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 scales associated with the shorter wavebands.  
+It is usually chosen to be 23~m through \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 the 
+particules concentration and it 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 
+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. 
+
+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 
+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 
+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 $z$-coordinate is preferred to $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. 
+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 wavebands chlorophyll-independant formulation (blue), 
+a chlorophyll-dependant 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}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1067,18 +1121,24 @@
 %--------------------------------------------------------------------------------------------------------------
 
-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+\rdt} &= T^{t-\rdt} + 2 \, \rdt  \ \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-\rdt} -2t^t+T^{t+\rdt}} \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{rn\_atfp} (\textbf{namelist} parameter). 
-Its default value is \np{rn\_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 
Index: /branches/DEV_r1826_DOC/DOC/TexFiles/Namelist/namsbc_rnf
===================================================================
--- /branches/DEV_r1826_DOC/DOC/TexFiles/Namelist/namsbc_rnf	(revision 2194)
+++ /branches/DEV_r1826_DOC/DOC/TexFiles/Namelist/namsbc_rnf	(revision 2195)
@@ -2,18 +2,20 @@
 &namsbc_rnf    !   runoffs namelist surface boundary condition
 !-----------------------------------------------------------------------
-!              !  file name  ! frequency (hours) ! variable   ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation !
-!              !             !  (if <0  months)  !   name     !  (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  !
-   sn_rnf      = 'runoff_1m' ,        -1.        , 'sorunoff' ,   .true.     , .true.  , 'yearly'  , ''       , ''
-   sn_cnf      = 'runoff_1m' ,         0.        , 'socoefr'  ,   .false.    , .true.  , 'yearly'  , ''       , ''
-   sn_sal_rnf  = 'runoff_1m' ,        -1.        , 'rosaline' ,   .true.     , .true.  , 'yearly'  , ''       , ''
-   sn_tmp_rnf  = 'runoff_1m' ,        -1.        , 'rotemper' ,   .true.     , .true.  , 'yearly'  , ''       , ''
-   sn_dep_rnf  = 'runoff_1m' ,         0.        , 'rodepth'  ,   .false.    , .true.  , 'yearly'  , ''       , ''
- 
-   cn_dir       = './'      !  root directory for the location of the runoff files
-   ln_rnf_emp   = .false.   !  runoffs included into precipitation field (T) or into a file (F)
-   ln_rnf_att   = .false.   !  runoffs attributes (temp, sal & depth) are specified in a file
-   ln_rnf_mouth = .false.   !  specific treatment at rivers mouths
-   rn_hrnf      =  15.e0    !  depth over which enhanced vertical mixing is used
+!               !   file name        ! frequency (h)   ! variable   ! time interp. !  clim  ! starting !
+!               !                    ! (if <0  months) !   name     !  (logical)   !  (0/1) !  record  !
+   sn_rnf     = 'runoff_1m_nomask.nc',      -12.       , 'sorunoff' ,   .true.     ,    1   ,    0     , ''   , ''
+   sn_cnf     = 'runoff_1m_nomask.nc',        0.       , 'socoefr'  ,   .false.    ,    1   ,    0     , ''   , ''
+   sn_s_rnf   = 'runoffs'            ,  24             , 'rosaline' ,   .true.     , .true. ,  'yearly', ''   , ''
+   sn_t_rnf   = 'runoffs'            ,  24             , 'rotemper' ,   .true.     , .true. ,  'yearly', ''   , ''
+   sn_dep_rnf = 'runoffs'            ,   0             , 'rodepth'  ,   .false.    , .true. ,  'yearly', ''   , ''
+!
+   cn_dir       = './'      ! directory in which the model is executed
+   ln_rnf_emp   =   .false. !  runoffs included into precipitation field (T) or into a file (F)
+   ln_rnf_mouth =   .false. !  specific treatment at rivers mouths
+   rn_hrnf      =   15.e0   !  depth over which enhanced vertical mixing is used
    rn_avt_rnf   =   1.e-3   !  value of the additional vertical mixing coef. [m2/s]
    rn_rfact     =   1.e0    !  multiplicative factor for runoff
+   ln_rnf_depth =  .false.  !  read in depth information for runoff
+   ln_rnf_temp  =  .false.  !  read in temperature information for runoff
+   ln_rnf_sal   =  .false.  !  read in salinity information for runoff
 /
