Context Navigation

← Previous Change
Next Change →

Chap_TRA.tex

Timestamp:

2008-02-09T15:13:48+01:00 (16 years ago)

Author:

Message:

trunk - update including Steven correction of the first 5 chapters (until DYN) and activation of Appendix A & B

File:

: 1 edited

trunk/DOC/BETA/Chapters/Chap_TRA.tex (modified) (35 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/DOC/BETA/Chapters/Chap_TRA.tex

-                      r781
+                      r817
 %        ==> 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
+\newpage
 $\ $\newline    % force a new ligne
 Using the representation described in Chap.~\ref{DOM}, several semi-discrete space
 forms of the tracer equations are available depending on the vertical
+Using the representation described in Chap.~\ref{DOM}, several semi-discrete
+space forms of the tracer equations are available depending on the vertical
 coordinate used and on the physics used. In all the equations presented
 here, the masking has been omitted for simplicity. One must be aware that
 …
 operator is used, the resulting field is multiplied by a mask.
+The two active tracers are potential temperature and salinity. Their prognostic equation can be summarized as follows:
+The two active tracers are potential temperature and salinity. Their prognostic
+equations can be summarized as follows:
 \begin{equation*}
 \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC}
 …
 \end{equation*}
+NXT stands for next, referring to the time-stepping. From left to right, the terms on the rhs of the tracer equations are the advection (ADV), the lateral diffusion (LDF), the vertical diffusion (ZDF), the contributions from the external forcings (SBC: Surface Boundary Condition, QSR: Solar Radiation penetration, and BBC: Bottom Boundary Condition), the contribution from the bottom boundary Layer (BBL) parametrisation, and an internal damping (DMP) term. The last four have been put inside brackets as they are optional. The external forcings and parameterizations require complex inputs and calculations (bulk formulae, estimation of mixing coefficients) that are carried out in modules of the SBC, LDF and ZDF categories and described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively. Note that \mdl{tranpc}, the non-penetrative convection module,  although (temporarily) located in the NEMO/OPA//TRA directory, is described with the model vertical physics (ZDF).
+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 CPP keys. For each equation term 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 found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
+The user has the option of extracting each tendency term on the rhs of the tracer equation (\key {trdtra} defined), as described in Chap.~\ref{MISC}.
+NXT stands for next, referring to the time-stepping. From left to right, the terms
+on the rhs of the tracer equations are the advection (ADV), the lateral diffusion
+(LDF), the vertical diffusion (ZDF), the contributions from the external forcings
+(SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC:
+Bottom Boundary Condition), the contribution from the bottom boundary Layer
+(BBL) parametrisation, and an internal damping (DMP) term. The terms QSR,
+BBC, BBL and DMP are optional. The external forcings and parameterizations
+require complex inputs and complex calculations (e.g. bulk formulae, estimation
+of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and
+described in chapters \S\ref{SBC}, \S\ref{LDF} and  \S\ref{ZDF}, respectively.
+Note that \mdl{tranpc}, the non-penetrative convection module,  although
+(temporarily) located in the NEMO/OPA/TRA directory, is described with the
+model vertical physics (ZDF).
+%%%
+\gmcomment{change the position of eosbn2 in the reference code}
+%%%
+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
+CPP keys. For each equation term 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
+found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
+The user has the option of extracting each tendency term on the rhs of the tracer
+equation for output (\key{trdtra} is defined), as described in Chap.~\ref{MISC}.
 % ================================================================
 % Tracer Advection
 % ================================================================
+\section{Tracer Advection (\mdl{traadv})}
+\section  [Tracer Advection (\textit{traadv})]
+      {Tracer Advection (\mdl{traadv})}
 \label{TRA_adv}
 %------------------------------------------nam_traadv-----------------------------------------------------
 …
 %-------------------------------------------------------------------------------------------------------------
 The advection tendency in flux form is the divergence of the advective
+The advection tendency of a tracer in flux form is the divergence of the advective
 fluxes. Its discrete expression is given by :
 \begin{equation} \label{Eq_tra_adv}
 …
 {w\;\tau _w } \right]
 \end{equation}
+which, in pure z-coordinate (\key{zco} defined), reduces to :
+where $\tau$ is either T or S. 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
 …
 e\nolimits_{3T} }\delta _k \left[ {w\;\tau _w } \right]
 \end{equation}
 as the vertical scale factors are function of $k$ only, and thus $e_{3u}
+since the vertical scale factors are functions of $k$ only, and thus $e_{3u}
 =e_{3v} =e_{3T} $.
 The flux form requires implicitly the use of the continuity equation:
+The flux form in \eqref{Eq_tra_adv} requires implicitly the use of the continuity equation:
 $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$
+using $\nabla \cdot \vect{U}=0)$ or $\partial _t e_3 +\nabla \cdot \vect{U}=0$ in variable volume case ($i.e.$ \key{vvl} defined). Therefore it is of
+(using $\nabla \cdot \vect{U}=0)$ or $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$
+ in variable volume case ($i.e.$ \key{vvl} defined). Therefore it is of
 paramount importance to design the discrete analogue of the advection
 tendency so that it is consistent with the continuity equation in order to
 enforce conservation properties of the continuous equations. In other words,
 by substituting $\tau$ by 1 in (\ref{Eq_tra_adv}) we recover discrete form of the
 continuity equation which is used to calculate the vertical velocity.
+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
+the continuity equation which is used to calculate the vertical velocity.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t] \label{Fig_adv_scheme}  \begin{center}
 \includegraphics[width=0.9\textwidth]{./Figures/Fig_adv_scheme.pdf}
+\caption{Schematic representation of some ways used to evaluate the tracer value at $u$-point and the amount of tracer exchanged between two neighbouring grid points. Upsteam biased scheme (ups): the upstream value is used and the black area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation is used and black + dark grey areas is exchanged. Monotonic upstream scheme for conservative laws (muscl):  a parabolic interpolation is used and black + dark grey + grey areas are exchanged. Second order scheme (cen2): the mean value is used and black + dark grey + grey + light grey areas are exchanged. Note that this illustration does not include the flux limiter used in ppm and muscl schemes.}
+\caption{Schematic representation of some ways used to evaluate the tracer value
+at $u$-point and the amount of tracer exchanged between two neighbouring grid
+points. Upsteam biased scheme (ups): the upstream value is used and the black
+area is exchanged. Piecewise parabolic method (ppm): a parabolic interpolation
+is used and the black and dark grey areas are exchanged. Monotonic upstream
+scheme for conservative laws (muscl):  a parabolic interpolation is used and black,
+dark grey and grey areas are exchanged. Second order scheme (cen2): the mean
+value is used and black, dark grey, grey and light grey areas are exchanged. Note
+that this illustration does not include the flux limiter used in ppm and muscl schemes.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+The advection schemes used in OPA differ by the choice made in space and
+time interpolation to define the value of the tracer at the velocity points (\ref{Fig_adv_scheme}).
+The key difference between the advection schemes used 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
 specified, since the normal velocity is zero there. At the sea surface the
+boundary condition depends on the type of sea surface chosen: (1) in
+rigid-lid formulation, $w=0$ at the surface, so the advective fluxes through the
+surface is zero ; (2) in non-linear free surface (variable volume case,
+\key{vvl} defined), convergence/divergence in the first ocean level moves
+up/down the free surface: there is no tracer advection through it so that
+the advective fluxes through the surface is also zero ; (3) in the 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 \mathord{\left/ {\vphantom {1
+}} \right. \kern-\nulldelimiterspace} 2} =T_{k=1} $. This boundary
+condition retains local conservation of tracer. Strict global conservation
+is not possible in linear free surface but is achieved to a good
+approximation since the non-conservative 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}).
+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 \S\ref{DYN}). Nevertheless, 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.
+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} $.
+\item [non-linear free surface:] (\key{vvl} is defined)
+convergence/divergence in the first ocean level moves the free surface
+up/down. There is no tracer advection through it so that the advective
+fluxes through the surface are also zero
+\end{description}
+In all cases, this boundary condition retains local conservation of tracer.
+Global conservation is obtained in both rigid-lid and non-linear free surface
+cases, but not in the linear free surface case. Nevertheless, in the latter
+case, it is achieved to a good approximation since the non-conservative
+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}).
+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
 setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The
 corresponding code can be found in \textit{traadv\_xxx.F90} module, where
 \textit{xxx} is a 3 or 4 letter acronym accounting for each scheme. Details
+corresponding code can be found in the \textit{traadv\_xxx.F90} module, where
+\textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. Details
 of the advection schemes are given below. The choice of an advection scheme
 is a complex matter which depends on the model physics, model resolution,
 type of tracer, as well as the issue of numerical cost.
+Note that (1) cen2, cen4 and TVD schemes require an explicit diffusion
+Note that
+(1) cen2, cen4 and TVD schemes require an explicit diffusion
 operator while the other schemes are diffusive enough so that they do not
+require additional diffusion ; (2) cen2, cen4, MUSCL2, and UBS are not
+\textit{positive} schemes, meaning false extrema are permitted. It is not recommended to use
+them on passive tracers ; (3) It is highly recommended to use the same
+advection-diffusion scheme on both active and passive tracers. In
+particular, if a source or sink of a passive tracer depends on a active one,
+the difference of treatment of active and passive tracers can create very
+nice-looking frontal structures that are pure numerical artefacts.
+require additional diffusion ;
+(2) cen2, cen4, MUSCL2, and UBS are not \textit{positive} schemes
+\footnote{negative values can appear in an initially strictly positive tracer field
+which is advected}
+, implying that false extrema are permitted. Their use is not recommended on passive tracers ;
+(3) It is highly recommended that the same advection-diffusion scheme is
+used on both active and passive tracers. Indeed, if a source or sink of a
+passive tracer depends on an active one, the difference of treatment of
+active and passive tracers can create very nice-looking frontal structures
+that are pure numerical artefacts.
 % -------------------------------------------------------------------------------------------------------------
 %        2nd order centred scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=T)}
+\subsection   [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})]
+         {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=.true.)}
 \label{TRA_adv_cen2}
 In the centred second order formulation, the tracer at velocity points is
 evaluated as the mean of the two neighbouring $T$-points. For example, in the
 $i$-direction :
+evaluated as the mean of the two neighbouring $T$-point values.
+For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_cen2}
 \tau _u^{cen2} =\overline T ^{i+1/2}
 \end{equation}
+The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$ but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to produce a sensible solution. The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
+Note that using 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.
+The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$
+but dispersive ($i.e.$ it may create false extrema). It is therefore notoriously
+noisy and must be used in conjunction with an explicit diffusion operator to
+produce a sensible solution. The associated time-stepping is performed using
+a leapfrog scheme in conjunction with an Asselin time-filter, so $T$ in
+(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
+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.
 % -------------------------------------------------------------------------------------------------------------
 %        4nd order centred scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4}=T)}
+\subsection   [$4^{nd}$ order centred scheme (cen4) (\np{ln\_traadv\_cen4})]
+         {$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 is evaluated
 at velocity points as the $4^{th}$ order interpolation of $T$, and thus use the
 four neighbouring $T$-points. For example, in the $i$-direction:
+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
+the four neighbouring $T$-points. For example, in the $i$-direction:
 \begin{equation} \label{Eq_tra_adv_cen4}
 \tau _u^{cen4}
 …
 Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme
 but a $4^{th}$ order evaluation of advective fluxes since the divergence of
 advective fluxes, (\ref{Eq_tra_adv}), is kept at $2^{nd}$ order. The ``$4^{th}$ order
 scheme'' denomination used in oceanographic literature is usually associated
+but a $4^{th}$ order evaluation of advective fluxes, since the divergence of
+advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order. The phrase ``$4^{th}$
+order scheme'' used in oceanographic literature is usually associated
 with the scheme presented here. Introducing a \textit{true} $4^{th}$ order advection
 scheme is feasible but, for consistency reasons, it requires changes in the
 discretisation of the tracer advection together with changes in both the
 continuity equation and the momentum advection.
+continuity equation and the momentum advection terms.
 A direct consequence of the pseudo-fourth order nature of the scheme is that
 it is not non-diffusive, i.e. the global variance of a tracer is not
 preserved through \textit{cen4}. Furthermore, it must be used in conjunction with an
+preserved using \textit{cen4}. Furthermore, it must be used in conjunction with an
 explicit diffusion operator to produce a sensible solution. The
 time-stepping is also performed using a leapfrog scheme in conjunction with
 an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
 At $T$-grid cell abutted to a boundary (coastline, bottom and surface), an
+At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an
 additional hypothesis must be made to evaluate $\tau _u^{cen4}$. This
 hypothesis usually reduces the order of the scheme. Here we choose to set
 the gradient of $T$ across the boundary to zero. Alternative conditions can be
 specified such as the reduction to a second order scheme for near boundary
 grid point.
+specified, such as a reduction to a second order scheme for these near boundary
+grid points.
 % -------------------------------------------------------------------------------------------------------------
 %        TVD scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=T)}
+\subsection   [Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]
+         {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=.true.)}
 \label{TRA_adv_tvd}
 In the Total Variance Dissipation (TVD) formulation, the tracer at velocity
 points is evaluated as a combination of upstream and centred scheme. For
+points is evaluated using a combination of an upstream and a centred scheme. For
 example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_tvd}
 …
 \end{equation}
 where $c_u$ is a flux limiter function taking values between 0 and 1. There
+exists many ways to define $c_u$., each correcponding to a different total
+variance decreasing scheme. The one chosen in OPA is described in \citet{Zalesak1979}. $c_u$ only departs from $1$ when the advective term produces a local
+extremum in the tracer field. The resulting scheme is quite 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.
+For stability reasons in (\ref{Eq_tra_adv_tvd}) $\tau _u^{cen2}$ is evaluated using the
+\textit{now} velocity (leap-frog environment: centred in time) while $\tau _u^{ups}$ is
+evaluated using the \textit{before} velocity (diffusive part: forward in time).
+exist many ways to define $c_u$, each correcponding to a different total
+variance decreasing scheme. The one chosen in \NEMO is described in \citet{Zalesak1979}. $c_u$ only departs from $1$ when the advective term
+produces a local extremum in the tracer field. The resulting scheme is quite
+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.
+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}$
+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
+used for the diffusive part.
 % -------------------------------------------------------------------------------------------------------------
 %        MUSCL scheme
 % -------------------------------------------------------------------------------------------------------------
 \subsection[MUSCL scheme  (\np{ln\_traadv\_muscl}=T)]
+\subsection[MUSCL scheme  (\np{ln\_traadv\_muscl})]
    {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_muscl}=T)}
 \label{TRA_adv_muscl}
 The Monotone Upstream Scheme for Conservative Laws (MUSCL) has been
 implemented by \citet{Levy2001}. 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 :
+implemented by \citet{Levy2001}. 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}
 …
    \end{aligned}    \right.
 \end{equation}
 where $\widetilde{\partial _i \tau}$ is the slope of the tracer
 on which a limitation is imposed to ensure the \textit{positive} character of the scheme.
 The time stepping is performed using a forward scheme, that is the \textit{before} tracer
 field is used to evaluate $\tau _u^{mus}$.
 For an ocean grid point abutted to land and where the ocean velocity is
 toward land, two choices are available: use of an upstream flux
 (\np{ln\_traadv\_muscl}=T) or use of second order flux
 (\np{ln\_traadv\_muscl2}=T). Note that the latter choice does not insure the
 \textit{positive} character of the scheme. Only the former can be used on both active and
 passive tracers.
+where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation
+is imposed to ensure the \textit{positive} character of the scheme.
+The time stepping is performed using a forward scheme, that is the \textit{before}
+tracer field is used to evaluate $\tau _u^{mus}$.
+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
+the \textit{positive} character of the scheme. Only the former can be used
+on both active and passive tracers.
 % -------------------------------------------------------------------------------------------------------------
 %        UBS scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Upstream Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=T)}
+\subsection   [Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs})]
+         {Upstream-Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=.true.)}
 \label{TRA_adv_ubs}
 The UBS advection scheme is an upstream biased third order scheme based on
 an upstream-biased parabolic interpolation. It is also known as Cell
+The UBS advection scheme is an upstream-biased third order scheme based on
+an upstream-biased parabolic interpolation. It is also known as the Cell
 Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective
 Kinematics). For example, in the $i$-direction :
 …
 advection scheme is similar to that reported in \cite{Farrow1995}.
 It is a relatively good compromise between accuracy and smoothness. It is
 not a \emph{positive} scheme meaning false extrema are permitted but the
+not a \emph{positive} scheme, meaning that false extrema are permitted, but the
 amplitude of such are significantly reduced over the centred second order
 method. Nevertheless it is not recommended to apply it to a passive tracer
 that requires positivity.
+method. Nevertheless it is not recommended that it should be applied to a passive
+tracer that requires positivity.
 The intrinsic diffusion of UBS makes its use risky in the vertical direction
 where the control of artificial diapycnal fluxes is of paramount importance.
+It has therefore been preferred to evaluate the vertical flux using the TVD
+scheme when \np{ln\_traadv\_ubs}=T.
+For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds to a
+second order centred scheme is evaluated using the \textit{now} velocity (centred in
+time) while the second term which is the diffusive part of the scheme, is
+evaluated using the \textit{before} velocity (forward in time. This is discussed by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK
+schemes only differ by one coefficient. Substituting 1/6 with 1/8 in
+(\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. This
+option is not available through a namelist parameter, since the 1/6
+coefficient is hard coded. Nevertheless it is quite easy to make the
+substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme
+NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can
+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)
+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 is discussed by \citet{Webb1998} in the context of the Quick
+advection scheme. UBS and QUICK
+schemes only differ by one coefficient. Replacing 1/6 with 1/8 in
+\eqref{Eq_tra_adv_ubs} leads to the QUICK advection scheme
+\citep{Webb1998}. This option is not available through a namelist
+parameter, since the 1/6 coefficient is hard coded. Nevertheless
+it is quite easy to make the substitution in the \mdl{traadv\_ubs} module
+and obtain a QUICK scheme.
+Note that :
+(1): When a high vertical resolution $O(1m)$ is used, the model stability can
 be controlled by vertical advection (not vertical diffusion which is usually
 solved using an implicit scheme). Computer time can be saved by using a
 time-splitting technique on vertical advection. This possibility have been
 implemented and validated in ORCA05-L301. It is not currently offered in the
+time-splitting technique on vertical advection. This case has been
+implemented and validated in ORCA05 with 301 levels. It is not available in the
 current reference version.
+NB 2 : In a forthcoming release four options will be proposed 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 \citet{Sacha2005} implementation of UBS in ROMS,
+or  \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an
+(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.
+following \citet{Sacha2005} implementation of UBS in ROMS, or  \textit{(d)}
+an UBS. The $3^{rd}$ case has dispersion properties similar to an
 eight-order accurate conventional scheme.
 NB 3 : It is straight forward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
+(3) : It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
 \begin{equation} \label{Eq_tra_adv_ubs2}
 \tau _u^{ubs} = \left\{  \begin{aligned}
 …
 - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i]
 \end{equation}
+\eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that the UBS scheme is based on the fourth order scheme to which is added an upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order part as stated above using \eqref{Eq_tra_adv_ubs} and also as it is coded in NEMO v2.3. Third, the diffusive term is in fact a biharmonic operator with a eddy coefficient with is simply proportional to the velocity: $A_u^{lm}= - \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note that the current version of NEMO uses \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_tra_adv_ubs2}.
+\eqref{Eq_tra_adv_ubs2} has several advantages. Firstly, it clearly reveals
+that the UBS scheme is based on the fourth order scheme to which an
+upstream-biased diffusion term is added. Secondly, this emphasises that the
+$4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has
+to be evaluated at the \emph{now} time step using \eqref{Eq_tra_adv_ubs}.
+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_tra_adv_ubs2}. This should be
+ changed in forthcoming release.
+ %%%
+ \gmcomment{the change in UBS scheme has to be done}
+ %%%
 % -------------------------------------------------------------------------------------------------------------
 %        QCK scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=T)}
+\subsection   [QUICKEST scheme (QCK) (\np{ln\_traadv\_qck})]
+         {QUICKEST scheme (QCK) (\np{ln\_traadv\_qck}=.true.)}
 \label{TRA_adv_qck}
 The Quadratic Upstream Interpolation for Convective Kinematics with
 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} is
 the third order Godunov scheme. It is associated with ULTIMATE QUICKEST
+Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979}
+is the third order Godunov scheme. It is associated with the ULTIMATE QUICKEST
 limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray
+(MERCATOR-ocean).
+The resulting scheme is quite expensive but \emph{positive}. It can be used on both active
+and passive tracers.
+(MERCATOR-ocean).
+The resulting scheme is quite expensive but \emph{positive}. It can be used on both active and passive tracers. Nevertheless, the intrinsic diffusion of QCK makes its use
+risky in the vertical direction where the control of artificial diapycnal fluxes is of
+paramount importance. Therefore the vertical flux is evaluated using the CEN2
+scheme. This no more ensure the positivity of the scheme. The use of TVD in the
+vertical direction as for the UBS case should be implemented to maintain the property.
 % -------------------------------------------------------------------------------------------------------------
 %        PPM scheme
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm}=T)}
+\subsection   [Piecewise Parabolic Method (PPM) (\np{ln\_traadv\_ppm})]
+         {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. As QCK scheme, it is associated with
+ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the current reference
+version.
+is based on a quadradic piecewise rebuilding. 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 2.3.
 % ================================================================
 % Tracer Lateral Diffusion
 % ================================================================
+\section{Tracer Lateral Diffusion (\mdl{traldf})}
+\section  [Tracer Lateral Diffusion (\textit{traldf})]
+      {Tracer Lateral Diffusion (\mdl{traldf})}
 \label{TRA_ldf}
 %-----------------------------------------nam_traldf------------------------------------------------------
 …
 %-------------------------------------------------------------------------------------------------------------
 The options available for lateral diffusion are laplacian (rotated or not)
 or biharmonic operators, that latter being more scale-selective (more
+The options available for lateral diffusion are a laplacian (rotated or not)
+or a biharmonic operator, the latter being more scale-selective (more
 diffusive at small scales). The specification of eddy diffusivity
 coefficients (either constant, variable in space and time) as well as the
 computation of the slope along which the operators act are performed in
+coefficients (either constant or variable in space and time) as well as the
+computation of the slope along which the operators act, are performed in the
 \mdl{ldftra} and \mdl{ldfslp} modules, respectively. This is described in Chap.~\ref{LDF}. The lateral diffusion of tracers is evaluated using a forward scheme,
 i.e. the tracers appearing in its expression are the \textit{before} tracers in time,
 except for the pure vertical component that appears when a tensor of
 rotation is used. This latter term is solved implicitly together with the
 vertical diffusion term (see \S\ref{DOM_nxt})
+$i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time,
+except for the pure vertical component that appears when a rotation tensor
+is used. This latter term is solved implicitly together with the
+vertical diffusion term (see \S\ref{DOM_nxt}).
 % -------------------------------------------------------------------------------------------------------------
 %        Iso-level laplacian operator
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Iso-level laplacian operator (\mdl{traldf\_lap}, \np{ln\_traldf\_lap}) }
+\subsection   [Iso-level laplacian operator (\textit{traldf\_lap} - \np{ln\_traldf\_lap})]
+         {Iso-level laplacian operator (\mdl{traldf\_lap} - \np{ln\_traldf\_lap}=.true.) }
 \label{TRA_ldf_lap}
+A laplacian diffusive operator (i.e. a harmonic operator) acting along the model surfaces is given by:
+A laplacian diffusion operator (i.e. a harmonic operator) acting along the model
+surfaces is given by:
 \begin{equation} \label{Eq_tra_ldf_lap}
 \begin{split}
 …
 \end{equation}
+This lateral operator is a \emph{horizontal} one (i.e. acting along geopotential surfaces) in
+$z$-coordinate with or without partial step, but it is simply an iso-level operator in $s$-coordinate.
+It is thus used when, in addition to \np{ln\_traldf\_lap}=T, we have \np{ln\_traldf\_level}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=F. In both cases, it significantly contributes to diapycnal mixing. It is therefore not recommended to use it.
+\textit{Notes} : In pure z-coordinate (\key{zco} defined), $e_{3u}=e_{3v}=e_{3T}$, so that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}).
+\textit{Notes} : In partial step $z$-coordinate (\np{ln\_zps}=T), tracers in horizontally adjacent cells are located at different depths in vicinity of the bottom. In this case, horizontal derivatives in (\ref{Eq_tra_ldf_lap}) at the bottom level require a specific treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
+This lateral operator is a \emph{horizontal} one ($i.e.$ acting along
+geopotential surfaces) in the $z$-coordinate with or without partial step,
+but is simply an iso-level operator in the $s$-coordinate.
+It is thus used when, in addition to \np{ln\_traldf\_lap}=.true., we have
+\np{ln\_traldf\_level}=.true., or both \np{ln\_traldf\_hor}=.true. and
+\np{ln\_zco}=.false.. In both cases, it significantly contributes to
+diapycnal mixing. It is therefore not recommended.
+Note that
+(1) In pure $z$-coordinate (\key{zco} is defined), $e_{3u}=e_{3v}=e_{3T}$, so
+that the vertical scale factors disappear from (\ref{Eq_tra_ldf_lap}).
+(2) In partial step $z$-coordinate (\np{ln\_zps}=.true.), tracers in horizontally
+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
+require a specific treatment. They are calculated in the \mdl{zpshde} module,
+described in \S\ref{TRA_zpshde}.
 % -------------------------------------------------------------------------------------------------------------
 %        Rotated laplacian operator
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Rotated laplacian operator (\mdl{traldf\_iso}, \np{ln\_traldf\_lap})}
+\subsection   [Rotated laplacian operator (\textit{traldf\_iso} - \np{ln\_traldf\_lap})]
+         {Rotated laplacian operator (\mdl{traldf\_iso} - \np{ln\_traldf\_lap}=.true.)}
 \label{TRA_ldf_iso}
+The general form of the second order lateral tracer subgrid scale physics (\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and $s$-coordinates:
+The general form of the second order lateral tracer subgrid scale physics
+(\ref{Eq_PE_zdf}) takes the following semi-discrete space form in $z$- and
+$s$-coordinates:
 \begin{equation} \label{Eq_tra_ldf_iso}
 …
  \end{split}
  \end{equation}
+where $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and the surface along which the diffusive operator acts ($i.e.$ horizontal or iso-neutral surfaces).
+ It is thus used when, in addition to \np{ln\_traldf\_lap}=T, we have \np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T. 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 using the mask technique (see \S\ref{LBC_coast}).
+The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical derivatives. For numerical stability, the vertical second derivative must be solved using the same implicit time scheme as those used in the vertical physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term is not computed in \mdl{traldf} module, but in \mdl{trazdf} module where, if iso-neutral mixing is used, the vertical mixing coefficient is simply increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.
+This formulation conserves the tracer but does not ensure the decrease of the tracer variance. Nevertheless the treatment performed on the slopes (see \S\ref{LDF}) allows to run safely without any additional background horizontal diffusion \citep{Guily2001}. An alternate scheme \citep{Griffies1998} which preserves both tracer and its variance is currently been tested in \NEMO.
+Note that in partial step $z$-coordinate (\np{ln\_zps}=T), the horizontal derivatives in \eqref{Eq_tra_ldf_iso} at the bottom level require a specific treatment. They are calculated in module zpshde, described in \S\ref{TRA_zpshde}.
+where $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
+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
+using the mask technique (see \S\ref{LBC_coast}).
+The operator in \eqref{Eq_tra_ldf_iso} involves both lateral and vertical
+derivatives. For numerical stability, the vertical second derivative must
+be solved using the same implicit time scheme as that used in the vertical
+physics (see \S\ref{TRA_zdf}). For computer efficiency reasons, this term
+is not computed in the \mdl{traldf} module, but in the \mdl{trazdf} module
+where, if iso-neutral mixing is used, the vertical mixing coefficient is simply
+increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.
+This formulation conserves the tracer but does not ensure the decrease
+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
+\citep{Griffies1998} which preserves both tracer and its variance is currently
+been tested in \NEMO.
+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}.
 % -------------------------------------------------------------------------------------------------------------
 %        Iso-level bilaplacian operator
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Iso-level bilaplacian operator (\mdl{traldf\_bilap}, \np{ln\_traldf\_bilap})}
+\subsection   [Iso-level bilaplacian operator (\textit{traldf\_bilap} - \np{ln\_traldf\_bilap})]
+         {Iso-level bilaplacian operator (\mdl{traldf\_bilap} - \np{ln\_traldf\_bilap}=.true.)}
 \label{TRA_ldf_bilap}
+The lateral fourth order operator formulation on tracers is obtained by applying (\ref{Eq_tra_ldf_lap}) twice. It requires an additional assumption on boundary conditions: first and third derivative terms normal to the coast are set to zero.
+It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have \np{ln\_traldf\_level}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=F. In both cases, it can contributes to diapycnal mixing even if it should be less than in the laplacian case. It is therefore not recommended to use it.
+\textit{Notes:} In the code, the bilaplacian routine does not call twice the laplacian
+routine but is rather a specific routine. This is due to the fact that we
+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,
+although less than in the laplacian case. It is therefore not recommended.
+Note that in the code, the bilaplacian routine does not call the laplacian
+routine twice but is rather a separate routine. This is due to the fact that we
 introduce the eddy diffusivity coefficient, A, in the operator as: $\nabla
 \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$ and instead of
+\cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$, instead of
 $-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$ where
 $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations
 …
 %        Rotated bilaplacian operator
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Rotated bilaplacian operator (\mdl{traldf\_bilapg}, \np{ln\_traldf\_bilap})}
+\subsection   [Rotated bilaplacian operator (\textit{traldf\_bilapg} - \np{ln\_traldf\_bilap})]
+         {Rotated bilaplacian operator (\mdl{traldf\_bilapg} - \np{ln\_traldf\_bilap}=.true.)}
 \label{TRA_ldf_bilapg}
+The lateral fourth order operator formulation on tracers is obtained by applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption on boundary conditions: first and third derivative terms normal to the coast, the bottom and the surface are set to zero.
+ It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have \np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T. Nevertheless, this rotated bilaplacian operator has never been seriously tested. No warranties that it is neither 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 unrotated operator.
+The lateral fourth order operator formulation on tracers is obtained by
+applying (\ref{Eq_tra_ldf_iso}) twice. It requires an additional assumption
+on boundary conditions: first and third derivative terms normal to the
+coast, the bottom and the surface are set to zero.
+It is used when, in addition to \np{ln\_traldf\_bilap}=T, we have
+\np{ln\_traldf\_iso}=T, or both \np{ln\_traldf\_hor}=T and \np{ln\_zco}=T.
+Nevertheless, this rotated bilaplacian operator has never been seriously
+tested. No warranties that it is neither 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
+unrotated operator.
 % ================================================================
 % Tracer Vertical Diffusion
 % ================================================================
+\section{Tracer Vertical Diffusion (\mdl{trazdf})}
+\section  [Tracer Vertical Diffusion (\textit{trazdf})]
+      {Tracer Vertical Diffusion (\mdl{trazdf})}
 \label{TRA_zdf}
 %--------------------------------------------namzdf---------------------------------------------------------
 …
 %--------------------------------------------------------------------------------------------------------------
+The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, based on a laplacian operator. The vertical diffusive operator given by (\ref{Eq_PE_zdf}) takes the following semi-discrete space form:
+The formulation of the vertical subgrid scale tracer physics is the same
+for all the vertical coordinates, and is based on a laplacian operator.
+The vertical diffusion operator given by (\ref{Eq_PE_zdf}) takes the
+following semi-discrete space form:
+(\ref{Eq_PE_zdf}) takes the following semi-discrete space form:
 \begin{equation} \label{Eq_tra_zdf}
 \begin{split}
 …
 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity
+coefficients on Temperature and Salinity, respectively. Generally, $A_w^{vT}=A_w^{vS}$ ecept when double diffusion mixing is parameterised (\key{zdfddm} defined). The way these coefficients can be evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when iso-neutral mixing is used, the both mixing coefficient are increased by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}.
+At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt must be specified. At the surface they are prescribed from the surface forcing (see \S\ref{TRA_sbc}), while at the bottom they are set to zero for heat and salt, unless a geothermal flux forcing is prescribed as a bottom boundary condition (\S\ref{TRA_bbc}).
+coefficients on Temperature and Salinity, respectively. Generally,
+$A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is
+parameterised (\key{zdfddm} is defined). The way these coefficients
+are evaluated is given in \S\ref{ZDF} (ZDF). Furthermore, when
+iso-neutral mixing is used, both mixing coefficients are increased
+by $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$
+to account for the vertical second derivative of \eqref{Eq_tra_ldf_iso}.
+At the surface and bottom boundaries, the turbulent fluxes of
+momentum, heat and salt must be specified. At the surface they
+are prescribed from the surface forcing (see \S\ref{TRA_sbc}),
+whilst at the bottom they are set to zero for heat and salt unless
+a geothermal flux forcing is prescribed as a bottom boundary
+condition (\S\ref{TRA_bbc}).
 The large eddy coefficient found in the mixed layer together with high
+vertical resolution implies a too restrictive constraint on the time step in
+explicit time stepping case (\np{ln\_zdfexp}=True). Therefore, the default implicit time stepping is generally preferred for the vertical diffusion as it overcomes the stability
+constraint. A forward time differencing scheme (\np{ln\_zdfexp}=T) 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 tracers and dynamics.
+vertical resolution implies that in the case of explicit time stepping
+(\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
+tracers and dynamics.
 % ================================================================
 …
 %        surface boundary condition
 % -------------------------------------------------------------------------------------------------------------
+\subsection{surface boundary condition (\mdl{trasbc})}
+\subsection   [Surface boundary condition (\textit{trasbc})]
+         {Surface boundary condition (\mdl{trasbc})}
 \label{TRA_sbc}
 …
 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}^C $, 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 d at the surface with the atmosphere, and $F_{wf}^C$ the forcing
+on the concentration associated with the water flux. The latter forcing has
+also two components: a direct effect of change in concentration associated
+with the tracer carried by the water flux, and an indirect concentration/dilution effect :
+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*}
+Two cases must be distinguished, the nonlinear free surface case (\key{vvl} 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, excepted 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 precipitations 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:
+\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}
 \begin{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 in 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 (excepted 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,
+%AMT  has this been corrected in the code?
+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 trasbc. Because of the hypothesis made for the temperature of precipitation and runoffs, for temperature $F_{wf}^e +F_{wf}^i =0$. The resulting forcing term for temperature is:
+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?
+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}
 …
 \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}).
+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
+EMP*S is not compensated by the advection of fluid through the top level: in
+the rigid lid case (contrary to the linear free surface), because \textit{w(k=1) = 0}. As a
+result, even if the budget of \textit{EMP} is zero in average over the whole ocean
+domain, the associated salt flux is not, as sea-surface salinity and \textit{EMP} are
+intrinsically correlated (high \textit{SSS} are found where evaporation is strong while
+low \textit{SSS} is usually associated with high precipitation or river runoff input).
+The $Q_{ns} $ and \textit{EMP} fields are defined and updated in \mdl{sbcmod} module (see
+\S\ref{SBC}).
+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}).
 % -------------------------------------------------------------------------------------------------------------
 %        Solar Radiation Penetration
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Solar Radiation Penetration (\mdl{traqsr})}
+\subsection   [Solar Radiation Penetration (\textit{traqsr})]
+         {Solar Radiation Penetration (\mdl{traqsr})}
 \label{TRA_qsr}
 %--------------------------------------------namqsr---------------------------------------------------------
+%--------------------------------------------namqsr--------------------------------------------------------
 \namdisplay{namqsr}
 %--------------------------------------------------------------------------------------------------------------
+When the penetrative solar radiation option is used (\np{ln\_flxqsr}=T, 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}=F). A term is thus added to the time evolution equation of temperature \eqref{Eq_PE_tra_T} while the surface boundary condition is modified to take into account only the non-penetrative part of the surface heat flux:
+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.).
+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
+modified to take into account only the non-penetrative part of the surface
+heat flux:
 \begin{equation} \label{Eq_PE_qsr}
 \begin{split}
 …
 \end{equation}
 A formulation including extinction coefficients is assumed for the downward irradiance $I$
 \citep{Paulson1977}:
+A formulation involving two extinction coefficients is assumed for the
+downward irradiance $I$ \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] classification: $\xi_1 = 0.35m$, $\xi_2 = 0.23m$ 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 modified in routine \mdl{traqsr}. Note that in $z$-coordinates, 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 central 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. At last, note that when the ocean is shallow (< 200~m), the part of the solar radiation can reach the ocean floor. In this case, we have chosen that all the radiation is absorbed at the last ocean level ($i.e.$ $I_w$ is masked).
+When coupling with a biology model (PISCES or LOBSTER), it is possible to calculate the light attenuation using information from the biology model. At the time of this writing, reading the light attenuation from a file is not implemented yet in the reference version.
+\colorbox{yellow}{case 4 bands and bio-coupling to add !!!}
+$\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.35m$, $\xi_2 = 0.23m$ 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
+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 !!!}  }
+%
 % -------------------------------------------------------------------------------------------------------------
 %        Bottom Boundary Condition
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Bottom Boundary Condition (\mdl{trabbc} + \key{bbc})}
+\subsection   [Bottom Boundary Condition (\textit{trabbc} - \key{bbc})]
+         {Bottom Boundary Condition (\mdl{trabbc} - \key{bbc})}
 \label{TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
 …
 \begin{figure}[!t] \label{Fig_geothermal}  \begin{center}
 \includegraphics[width=1.0\textwidth]{./Figures/Fig_TRA_geoth.pdf}
+\caption{Geothermal Heat flux (in $mW.m^{-2}$) as inferred from the age of the sea floor and the formulae of \citet{Stein1992}.}
+\caption{Geothermal Heat flux (in $mW.m^{-2}$) as inferred from the age
+of the sea floor and the formulae of \citet{Stein1992}.}
 \end{center}   \end{figure}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 Usually it is considered that there is no exchange of heat nor salt through
 the ocean bottom, i.e. a no flux boundary condition is applied on active
 tracers at the bottom. This is the default option in NEMO, and it is
 implemented using the masking technique. Nevertheless, there exists a
 non-zero heat flux across the seafloor that is associated with the solid
 earth cooling. This flux is weak compared with surface fluxes --- a mean
 global value of $\sim0.1\;W/m^2$ \citep{Stein1992} --- but it is
 systematically positive and it acts only on the densest water masses. Taking
 this flux into account in a global ocean model increases by a few Sverdrups
+Usually it is assumed that there is no exchange of heat or salt through
+the ocean bottom, $i.e.$ a no flux boundary condition is applied on active
+tracers at the bottom. This is the default option in \NEMO, and it is
+implemented using the masking technique. Hoever, there is a
+non-zero heat flux across the seafloor that is associated with solid
+earth cooling. This flux is weak compared to surface fluxes (a mean
+global value of $\sim0.1\;W/m^2$ \citep{Stein1992}), but it is
+systematically positive and acts on the densest water masses. Taking
+this flux into account in a global ocean model increases
 the deepest overturning cell (i.e. the one associated with the Antarctic
 Bottom Water).
 The presence or not of a geothermal heating is controlled by the namelist
 parameter  \np{ngeo\_flux}. Set to 1, a constant geothermal heatingis
 introducted which value is given by the \np{ngeo\_flux\_const}, also a
 namelist parameter. Set to 2, a spatially varying geothermal heat flux is
 introducted which is provided in the geothermal\_heating.nc NetCDF
 file (Fig.\ref{Fig_geothermal}).
+Bottom Water) by a few Sverdrups.
+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
+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,
+a spatially varying geothermal heat flux is introduced which is provided
+in the geothermal\_heating.nc NetCDF file (Fig.\ref{Fig_geothermal}).
 % ================================================================
 % Bottom Boundary Layer
 % ================================================================
+\section{Bottom Boundary Layer (\mdl{trabbl} + \key{bbl\_diff} or \key{bbl\_adv})}
+\section  [Bottom Boundary Layer (\textit{trabbl}, \textit{trabbl\_adv} )]
+      {Bottom Boundary Layer (\mdl{trabbl}, \mdl{trabbl\_adv})}
 \label{TRA_bbl}
 %--------------------------------------------nambbl---------------------------------------------------------
 …
 %--------------------------------------------------------------------------------------------------------------
 In z-coordinate configuration, the bottom topography is represented as a
+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 those gravity plumes is critical to determine the
 density and volume flux of the densest waters of the ocean, such as the
 Antarctic Bottom water, or the North Atlantic Deep Water. $z$-coordinate
 models tend to overestimate the entrainment because the gravity flow is
+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 $s$-coordinate when
+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 parameterization first introduced by
+a sill \citep{Willebrand2001}, and the thickness of the plume is not resolved.
+The idea of the bottom boundary layer (BBL) parameterization first introduced by
 \citet{BeckDos1998} is to allow a direct communication between
 two adjacent bottom cells at varying level, whenever the densest water is
 located above the less dense water. The communication can be by diffusive
 fluxes (diffusive BBL), advective fluxes (advective BBL) or both. Only
 tracers are modified, not the velocities. Implementing a BBL
 parameterization for momentum is a more complex problem because of the
+pressure gradient errors.
+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
+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}.
 % -------------------------------------------------------------------------------------------------------------
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Diffusive Bottom Boundary layer (\mdl{trabbl})}
+\subsection{Diffusive Bottom Boundary layer (\key{bbl\_diff})}
 \label{TRA_bbl_diff}
+The lateral diffusivity $A_l^\sigma $ in the BBL can be prescribed with a
+spatial dependence, e.g., in the conditional form
+\begin{equation} \label{Eq_tra_bbl}
+When applying sigma-diffusion (\key{trabbl} is defined), the diffusive flux between
+two adjacent cells living at the ocean bottom 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
+\begin{equation} \label{Eq_tra_bbl_coef}
 A_l^\sigma (i,j,t)=\left\{ {\begin{array}{l}
  \mbox{large}\quad if\;\nabla \rho \cdot \nabla H<0 \\
+ A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\
  \\
 \quad \quad \;\,\mbox{otherwise} \\
  \end{array}} \right.
 \end{equation}
+The large value of the coefficient when the diffusive BBL is active is given
+by the namelist parameter \np{atrbbl.}
+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.
 % -------------------------------------------------------------------------------------------------------------
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Advective Bottom Boundary Layer (\mdl{trabb\_adv})}
+\subsection   {Advective Bottom Boundary Layer (\key{bbl\_adv})}
 \label{TRA_bbl_adv}
+Implemented in NEMO v2.
+\colorbox{yellow} {Documentation to be added here }
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t] \label{Fig_bbl}  \begin{center}
+\includegraphics[width=1.0\textwidth]{./Figures/Fig_BBL_adv.pdf}
+\caption{Advective Bottom Boundary Layer.}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+%%%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 \\
+ \\
+\quad \quad \;\,\mbox{otherwise} \\
+ \end{array}} \right.
+\end{equation}
 % ================================================================
 % Tracer damping
 % ================================================================
+\section{Tracer damping (\mdl{tradmp})}
+\section  [Tracer damping (\textit{tradmp})]
+      {Tracer damping (\mdl{tradmp})}
 \label{TRA_dmp}
 %--------------------------------------------namdmp-----------------------------------------------------
 …
 %--------------------------------------------------------------------------------------------------------------
 In some applications it can be useful to add a Newtonian damping term in the
 temperature and salinity equations:
+In some applications it can be useful to add a Newtonian damping term
+into the temperature and salinity equations:
 \begin{equation} \label{Eq_tra_dmp}
 \begin{split}
 …
  \end{split}
  \end{equation}
+where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ 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 ($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
+\rou{dtacof} routine also located in  \mdl{tradmp}.
+The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)} the specification of
+the boundary conditions along artificial walls of a limited domain basin and
+\textit{(b)} the computation of the velocity field associated with a given $T$-$S$ field
+(for example to build the initial state of a prognostic simulation, or to
+use the resulting velocity 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 a few day time scale) 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 to find the velocity
+field consistent with the model dynamics while having a $T$-$S$ field close to a
+given climatology field ($T_o -S_o$). The time scale associated with
+$S_o$ is generally not a constant but spatially varying in order to respect
+some considerations. 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} as those two regions have a small time scale of adjustment,
+while smaller $S_o$ are used in the deep ocean where the typical time scale
+is long \citep{Sarmiento1982}. In addition it is reduced (and even
+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 \np{nmldmp}
+(\textbf{namelist} \np{nmldmp}parameter).
+The robust diagnostic method is very efficient to prevent the temperature drift in intermediate waters but it produces artificial sources of heat and salt within the ocean. It has also undesirable effects on the ocean convection. It tends to prevent deep convection and subsequent deep-water formation by stabilising too much the water columns.
+An example of computation of $S_o$ for robust diagnostic experiments with the ORCA2 model is provided in the \mdl{tradmp} module (subroutines \rou{dtacof} and \rou{cofdis} which compute coefficient and the distance to the bathymetry, respectively). Those routines are provided as examples and can be customised by the user.
+where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$
+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
+($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
+\rou{dtacof} also located in module \mdl{tradmp}.
+The two main cases in which \eqref{Eq_tra_dmp} is used are \textit{(a)}
+the specification of the boundary conditions along artificial walls of a
+limited domain basin and \textit{(b)} the computation of the velocity
+field associated with a given $T$-$S$ field (for example to build the
+initial state of a prognostic simulation, or to use the resulting velocity
+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
+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
+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
+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}.
+The robust diagnostic method is very efficient in preventing temperature
+drift in intermediate waters but it produces artificial sources of heat and salt
+within the ocean. It also has undesirable effects on the ocean convection.
+It tends to prevent deep convection and subsequent deep-water formation,
+by stabilising the water column too much.
+An example of the computation of $S_o$ for robust diagnostic experiments
+with the ORCA2 model is provided in the \mdl{tradmp} module
+(subroutines \rou{dtacof} and \rou{cofdis} which compute the coefficient
+and the distance to the bathymetry, respectively). These routines are
+provided as examples and can be customised by the user.
 % ================================================================
 % Tracer time evolution
 % ================================================================
+\section{Tracer time evolution (\mdl{tranxt})}
+\section  [Tracer time evolution (\textit{tranxt})]
+      {Tracer time evolution (\mdl{tranxt})}
 \label{TRA_nxt}
 %--------------------------------------------namdom-----------------------------------------------------
 …
 %--------------------------------------------------------------------------------------------------------------
+The general framework of dynamics 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 leap-frog scheme,
+$i.e.$ a three level centred time scheme associated with a Asselin time
+filter (cf. \S\ref{DOM_nxt}):
 \begin{equation} \label{Eq_tra_nxt}
 \begin{split}
 …
 \end{equation}
+where $\text{RHS}_T$ is the right hand side of the temperature equation, the subscript $f$ denotes
+filtered values and $\gamma$ is the Asselin coefficient. $\gamma$ is initialized as \np{atfp} (\textbf{namelist} parameter). Its default value is \np{atfp=0.1}.
+When the vertical mixing is solved implicitly, the update of the next tracer
+fields is done in module \mdl{trazdf}. In that case only the swap of arrays
+and the Asselin filtering is done in \mdl{tranxt} module.
+In order to prepare the computation of the next time step, a swap of tracer arrays is performed: $T^{t-\Delta t} = T^t$ and $T^t = T_f$.
+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}.
+When the vertical mixing is solved implicitly, the update of the \textit{next} tracer
+fields is done in module \mdl{trazdf}. In this case only the swapping of arrays
+and the Asselin filtering is done in the \mdl{tranxt} module.
+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$.
 % ================================================================
 % Equation of State (eosbn2)
 % ================================================================
+\section{Equation of State (\mdl{eosbn2}) }
+\section  [Equation of State (\textit{eosbn2}) ]
+      {Equation of State (\mdl{eosbn2}) }
 \label{TRA_eosbn2}
 %--------------------------------------------nameos-----------------------------------------------------
 …
 \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 is the equation of state 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. Even though this computation is fully vectorised, it is quite time consuming ($15$ to $20${\%} of the total CPU time) as 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 chosen the \citet{JackMcD1995} equation of state for seawater. It allows the computation of the \textit{in situ} ocean density directly as a function of \textit{potential} temperature relative to the sea surface (an OPA variable), the practical salinity (another OPA 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 the same expression except that the values of the various coefficients have been adjusted by \citet{JackMcD1995} in order to use directly 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 over $\rho_o$, a reference volumic mass (\textit{rau0} defined in \mdl{phycst}, usually $rau0= 1,020~Kg/m^3$). The default option (\np{neos}=0) is the \citet{JackMcD1995} equation of state. It is highly recommended to use it. Nevertheless, 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):
+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):
 \begin{equation} \label{Eq_tra_eos_linear}
 \begin{aligned}
 …
 \end{aligned}
 \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 \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 used as such.
+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
+used as such.
 % -------------------------------------------------------------------------------------------------------------
 …
 \label{TRA_bn2}
+An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a} frequency) is of paramount importance as it is used in several ocean parameterisations (namely TKE, KPP, Richardson number dependent vertical diffusion, enhanced vertical diffusion, non-penetrative convection, iso-neutral diffusion). In particular, one must be aware that $N^2$ has to be computed with an \textit{in situ} reference. The expression of $N^2$ depends on the type of equation of state used (\np{neos} namelist parameter).
+An accurate computation of the ocean stability (i.e. of $N$, the brunt-Vais\"{a}l\"{a}
+ frequency) is of paramount importance as it is used in several ocean
+ parameterisations (namely TKE, KPP, Richardson number dependent
+ vertical diffusion, enhanced vertical diffusion, non-penetrative convection,
+ 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}
+polynomial expression is used with the pressure in decibar approximated by
+the depth in meters:
+\begin{multline} \label{Eq_tra_bn2}
+N^2 = \frac{g}{e_{3w}} \; \beta ( \overline{T}^{\,k+1/2},\widetilde{S},z_w )   \\
+      \left\{  \alpha / \beta ( \overline{T}^{\,k+1/2},\widetilde{S},z_w )
+         \ \delta_{k+1/2}[T]     - \delta_{k+1/2}[S]
+       \right\}
+\end{multline}
+where $T$ is the \textit{potential} temperature, $\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$ a salinity anomaly, and $\alpha$ ($\beta\,$) the thermal (haline) expansion coefficient. Both $\alpha$ and $\beta$ depend on \textit{potential} temperature, salinity which are averaged at $w$-points prior to the computation.
+When a linear equation of state is used (\np{neos}=1 or 2, \eqref{Eq_tra_bn2} reduces to:
+polynomial expression is used (with the pressure in decibar approximated by
+the depth in meters):
+\begin{equation} \label{Eq_tra_bn2}
+N^2 = \frac{g}{e_{3w}} \; \beta   \
+      \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.
+Note that both $\alpha$ and $\beta$ depend on \textit{potential}
+temperature and salinity which are averaged at $w$-points prior
+to the computation instead of being computed at $T$-points and
+then averaged to $w$-points.
+When a linear equation of state is used (\np{neos}=1 or 2,
+\eqref{Eq_tra_bn2} reduces to:
 \begin{equation} \label{Eq_tra_bn2_linear}
 N^2 = \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right)
 \end{equation}
 where $\alpha$ and $\beta $ are the constant coefficients used to defined the linear
 equation of state \eqref{Eq_tra_eos_linear}.
+where $\alpha$ and $\beta $ are the constant coefficients used to
+defined the linear equation of state \eqref{Eq_tra_eos_linear}.
 % -------------------------------------------------------------------------------------------------------------
 %        Specific Heat
 % -------------------------------------------------------------------------------------------------------------
+\subsection{Specific Heat (\textit{rcp}, \mdl{phycst})}
+\subsection   [Specific Heat (\textit{phycst})]
+         {Specific Heat (\mdl{phycst})}
 \label{TRA_adv_ldf}
+The specific heat of sea water, $C_p$, is a function of temperature, salinity and pressure \citep{UNESCO1983}. It is only used in the model to convert surface heat fluxes into surface temperature increase, thus the pressure dependence is neglected. The dependence on $T$ and $S$ is weak. For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$ when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$. Its computer name is \textit{rcp} and its value is set in \mdl{phycst} module.
+The specific heat of sea water, $C_p$, is a function of temperature, salinity
+and pressure \citep{UNESCO1983}. It is only used in the model to convert
+surface heat fluxes into surface temperature increase and so the pressure
+dependence is neglected. The dependence on $T$ and $S$ is weak.
+For example, with $S=35~psu$, $C_p$ increases from $3989$ to $4002$
+when $T$ varies from -2~\degres C to 31~\degres C. Therefore, $C_p$ has
+been chosen as a constant: $C_p=4.10^3~J\,Kg^{-1}\,\degres K^{-1}$.
+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 (\mdl{ocfzpt})}
+\subsection   [Freezing Point of Seawater (\textit{ocfzpt})]
+         {Freezing Point of Seawater (\mdl{ocfzpt})}
 \label{TRA_fzp}
 …
 \end{equation}
+\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) has 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 0 or 1 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.
+\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}
+%%%
 % ================================================================
 % Horizontal Derivative in zps-coordinate
 % ================================================================
+\section{Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})}
+\section  [Horizontal Derivative in \textit{zps}-coordinate (\textit{zpshde})]
+      {Horizontal Derivative in \textit{zps}-coordinate (\mdl{zpshde})}
 \label{TRA_zpshde}
+With partial bottom cells (\np{ln\_zps}=T), tracers in horizontally adjacent cells generally live at different depths. Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure gradient (\mdl{dynhpg} module). Before taking horizontal gradients between the tracers next to the bottom, a linear interpolation is used to approximate the deeper tracer as if it actually lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}). For example on temperature in the i-direction, the needed interpolated temperature, $\widetilde{T}$, is:
+\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
+adjacent cells live at different depths. Horizontal gradients of tracers are needed
+for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure
+gradient (\mdl{dynhpg} module) to be active.
+\gmcomment{STEVEN from gm : question: not sure of  what -to be active- means}
+Before taking horizontal gradients between the tracers next to the bottom, a linear
+interpolation in the vertical is used to approximate the deeper tracer as if it actually
+lived at the depth of the shallower tracer point (Fig.~\ref{Fig_Partial_step_scheme}).
+For example, for temperature in the $i$-direction the needed interpolated
+temperature, $\widetilde{T}$, is:
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!p] \label{Fig_Partial_step_scheme}  \begin{center}
 \includegraphics[width=0.9\textwidth]{./Figures/Partial_step_scheme.pdf}
 \caption{ Discretisation of horizontal derivative and mean of tracers in z-partial step coordinate (\np{ln\_zps}=T) 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 derivative is then given by: $\delta _{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and the mean 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}
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 …
             \end{aligned}   \right.
 \end{equation*}
+and the resulting formulation of horizontal derivative and horizontal mean value of $T$ at $U$-point are:
+and the resulting forms for the horizontal difference and the horizontal average
+value of $T$ at a $U$-point are:
 \begin{equation} \label{Eq_zps_hde}
 \begin{aligned}
 …
 \end{equation}
+The computation of horizontal derivative of tracers as well as of density is performed once for all at each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. It has to be emphasized that the computation of the interpolated density, $\widetilde{\rho}$, is not identical to the one of $T$ and $S$. Instead of forming a linear approximation of density, we compute $\widetilde{\rho }$ from the interpolated value of $T$ and $S$, and the pressure of at $u$-point (in the equation of state pressure is approximated by depth, see \S\ref{TRA_eos} ) :
+The computation of horizontal derivative of tracers as well as of density is
+performed once for all at each time step in \mdl{zpshde} module and stored
+in shared arrays to be used when needed. It has to be emphasized that the
+procedure used to compute the interpolated density, $\widetilde{\rho}$, is not
+the same as that used for $T$ and $S$. Instead of forming a linear approximation
+of density, we compute $\widetilde{\rho }$ from the interpolated values of $T$
+and $S$, and the pressure at a $u$-point (in the equation of state pressure is
+approximated by depth, see \S\ref{TRA_eos} ) :
 \begin{equation} \label{Eq_zps_hde_rho}
 \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u })
 …
 \end{equation}
+This is a much better approximation as the variation of $\rho$ with depth (and thus pressure) is highly non-linear with a true equation of state and thus is badly approximated with a linear interpolation. This approximation is used to compute both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral surfaces (\S\ref{LDF_slp})
+\textit{Notes}: in almost all the advection schemes presented in this Chapter, both mean and derivative operators appear. Yet, it has been chosen not to use \eqref{Eq_zps_hde} in those schemes.: contrary to diffusion and pressure gradient computation, no correction for partial steps is applied for advection.The main motivation was to preserve the domain averaged mean variance of the field advected when using $2^{nd}$ order centred scheme. Sensitivity of the advection schemes to the way horizontal means are performed in the vicinity of partial cells should be further investigated in a near future.
+This is a much better approximation as the variation of $\rho$ with depth (and
+thus pressure) is highly non-linear with a true equation of state and thus is badly
+approximated with a linear interpolation. This approximation is used to compute
+both the horizontal pressure gradient (\S\ref{DYN_hpg}) and the slopes of neutral
+surfaces (\S\ref{LDF_slp})
+Note that in almost all the advection schemes presented in this Chapter, both
+averaging and differencing operators appear. Yet \eqref{Eq_zps_hde} has not
+been used in these schemes: in contrast to diffusion and pressure gradient
+computations, no correction for partial steps is applied for advection. The main
+motivation is to preserve the domain averaged mean variance of the advected
+field when using the $2^{nd}$ order centred scheme. Sensitivity of the advection
+schemes to the way horizontal averages are performed in the vicinity of partial
+cells should be further investigated in the near future.
+%%%
+\gmcomment{gm :   this last remark has to be done}
+%%%

Note: See TracChangeset for help on using the changeset viewer.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

Changeset 817 for trunk/DOC/BETA/Chapters/Chap_TRA.tex

Legend:

trunk/DOC/BETA/Chapters/Chap_TRA.tex

Download in other formats: