Context Navigation

← Previous Change
Next Change →

Chap_TRA.tex

Timestamp:

2016-11-28T17:04:10+01:00 (7 years ago)

Author:

emanuelaclementi

Message:

ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161

File:

: 1 edited

branches/2016/dev_INGV_UKMO_2016/DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (34 diffs)

Legend:

: Unmodified
: Added
: Removed

branches/2016/dev_INGV_UKMO_2016/DOC/TexFiles/Chapters/Chap_TRA.tex

-                      r5102
+                      r7351
+% ================================================================
+% Chapter 1 � Ocean Tracers (TRA)
+\documentclass[NEMO_book]{subfiles}
+\begin{document}
+% ================================================================
+% Chapter 1 ——— Ocean Tracers (TRA)
 % ================================================================
 \chapter{Ocean Tracers (TRA)}
 …
 (BBL) parametrisation, and an internal damping (DMP) term. The terms QSR,
 BBC, BBL and DMP are optional. The external forcings and parameterisations
 require complex inputs and complex calculations (e.g. bulk formulae, estimation
+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}
+%%%
+Note that \mdl{tranpc}, the non-penetrative convection module, although
+located in the NEMO/OPA/TRA directory as it directly modifies the tracer fields,
+is described with the model vertical physics (ZDF) together with other available
+parameterization of convection.
 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
+the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and
 freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}).
 The different options available to the user are managed by namelist logicals or
 CPP keys. For each equation term \textit{ttt}, the namelist logicals are \textit{ln\_trattt\_xxx},
+The different options available to the user are managed by namelist logicals or CPP keys.
+For each equation term  \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx},
 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme.
 The CPP key (when it exists) is \textbf{key\_trattt}. The equivalent code can be
 found in the \textit{trattt} or \textit{trattt\_xxx} module, in the NEMO/OPA/TRA directory.
 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}.
+The CPP key (when it exists) is \textbf{key\_traTTT}. The equivalent code can be
+found in the \textit{traTTT} or \textit{traTTT\_xxx} module, in the NEMO/OPA/TRA directory.
+The user has the option of extracting each tendency term on the RHS of the tracer
+equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{DIA}.
 $\ $\newline    % force a new ligne
 …
 %-------------------------------------------------------------------------------------------------------------
+The advection tendency of a tracer in flux form is the divergence of the advective
+fluxes. Its discrete expression is given by :
+When considered ($i.e.$ when \np{ln\_traadv\_NONE} is not set to \textit{true}),
+the advection tendency of a tracer is expressed in flux form,
+$i.e.$ as the divergence of the advective fluxes. Its discrete expression is given by :
 \begin{equation} \label{Eq_tra_adv}
 ADV_\tau =-\frac{1}{b_t} \left(
 …
 implicitly requires the use of the continuity equation. Indeed, it is obtained
 by using the following equality : $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$
 which results from the use of the continuity equation, $\nabla \cdot \vect{U}=0$ or
 $ \partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ in constant volume or variable volume case, respectively.
+which results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$
+(which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, $i.e.$ \np{ln\_linssh}=true).
 Therefore it is of paramount importance to design the discrete analogue of the
 advection tendency so that it is consistent with the continuity equation in order to
 enforce the conservation properties of the continuous equations. In other words,
 by replacing $\tau$ by the number 1 in (\ref{Eq_tra_adv}) we recover the discrete form of
+by setting $\tau = 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]    \begin{center}
 \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Fig_adv_scheme.pdf}
+\includegraphics[width=0.9\textwidth]{Fig_adv_scheme}
 \caption{   \label{Fig_adv_scheme}
 Schematic representation of some ways used to evaluate the tracer value
 …
 boundary condition depends on the type of sea surface chosen:
 \begin{description}
 \item [linear free surface:] the first level thickness is constant in time:
+\item [linear free surface:] (\np{ln\_linssh}=true) 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
 …
 $\left. {\tau _w } \right|_{k=1/2} =T_{k=1} $, $i.e.$
 the product of surface velocity (at $z=0$) by the first level tracer value.
 \item [non-linear free surface:] (\key{vvl} is defined)
+\item [non-linear free surface:] (\np{ln\_linssh}=false)
 convergence/divergence in the first ocean level moves the free surface
 up/down. There is no tracer advection through it so that the advective
 …
 \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
+Global conservation is obtained in non-linear free surface case,
+but \textit{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{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}).
+height, two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}.
 The velocity field that appears in (\ref{Eq_tra_adv}) and (\ref{Eq_tra_adv_zco})
+is the centred (\textit{now}) \textit{eulerian} ocean velocity (see Chap.~\ref{DYN}).
+When eddy induced velocity (\textit{eiv}) parameterisation is used it is the \textit{now}
+\textit{effective} velocity ($i.e.$ the sum of the eulerian and eiv velocities) which is used.
+The choice of an advection scheme is made in the \textit{\ngn{nam\_traadv}} namelist, by
+setting to \textit{true} one and only one of the logicals \textit{ln\_traadv\_xxx}. The
+corresponding code can be found in the \textit{traadv\_xxx.F90} module, where
+\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 the centred (\textit{now}) \textit{effective} ocean velocity, $i.e.$ the \textit{eulerian} velocity
+(see Chap.~\ref{DYN}) plus the eddy induced velocity (\textit{eiv})
+and/or the mixed layer eddy induced velocity (\textit{eiv})
+when those parameterisations are used (see Chap.~\ref{LDF}).
+Several tracer advection scheme are proposed, namely
+a $2^{nd}$ or $4^{th}$ order centred schemes (CEN),
+a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT),
+a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL),
+a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), and
+a Quadratic Upstream Interpolation for Convective Kinematics with
+Estimated Streaming Terms scheme (QUICKEST).
+The choice is made in the \textit{\ngn{namtra\_adv}} namelist, by
+setting to \textit{true} one of the logicals \textit{ln\_traadv\_xxx}.
+The 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.
+By default ($i.e.$ in the reference namelist, \ngn{namelist\_ref}), all the logicals
+are set to \textit{false}. If the user does not select an advection scheme
+in the configuration namelist (\ngn{namelist\_cfg}), the tracers will \textit{not} be advected !
+Details of the advection schemes are given below. The choosing 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
+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
+type of tracer, as well as the issue of numerical cost. In particular, we note that
+(1) CEN and FCT schemes require an explicit diffusion operator
+while the other schemes are diffusive enough so that they do not necessarily need additional diffusion ;
+(2) CEN and UBS are not \textit{positive} schemes
 \footnote{negative values can appear in an initially strictly positive tracer field
 which is advected}
 …
 % -------------------------------------------------------------------------------------------------------------
+%        2nd and 4th order centred schemes
+% -------------------------------------------------------------------------------------------------------------
+\subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})]
+            {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
+\label{TRA_adv_cen}
 %        2nd order centred scheme
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2})]
+         {$2^{nd}$ order centred scheme (cen2) (\np{ln\_traadv\_cen2}=true)}
+\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$-point values.
+The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~=~\textit{true}.
+Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
+and vertical direction by setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$.
+CEN implementation can be found in the \mdl{traadv\_cen} module.
+In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points
+is evaluated as the mean of the two neighbouring $T$-point values.
 For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_cen2}
 …
 \end{equation}
 The scheme is non diffusive ($i.e.$ it conserves the tracer variance, $\tau^2)$
+CEN2 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. The centered second
+order advection is computed in the \mdl{traadv\_cen2} module. In this module,
+it is advantageous to combine the \textit{cen2} scheme with an upstream scheme
+in specific areas which require a strong diffusion in order to avoid the generation
+of false extrema. These areas are the vicinity of large river mouths, some straits
+with coarse resolution, and the vicinity of ice cover area ($i.e.$ when the ocean
+temperature is close to the freezing point).
+This combined scheme has been included for specific grid points in the ORCA2
+and ORCA4 configurations only. This is an obsolescent feature as the recommended
+advection scheme for the ORCA configuration is TVD (see  \S\ref{TRA_adv_tvd}).
+Note that using the cen2 scheme, the overall tracer advection is of second
+(\ref{Eq_tra_adv_cen2}) is the \textit{now} tracer value.
+Note that using the CEN2, 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. \gmcomment{Note also that ... blah, blah}
+% -------------------------------------------------------------------------------------------------------------
+have this order of accuracy.
 %        4nd order centred scheme
+% -------------------------------------------------------------------------------------------------------------
+\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 values are
+evaluated at velocity points as a $4^{th}$ order interpolation, and thus depend on
+the four neighbouring $T$-points. For example, in the $i$-direction:
+In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as
+a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points.
+For example, in the $i$-direction:
 \begin{equation} \label{Eq_tra_adv_cen4}
 \tau _u^{cen4}
 =\overline{   T - \frac{1}{6}\,\delta _i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2}
 \end{equation}
+Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme
+In the vertical direction (\np{nn\_cen\_v}=$4$), a $4^{th}$ COMPACT interpolation
+has been prefered \citep{Demange_PhD2014}.
+In the COMPACT scheme, both the field and its derivative are interpolated,
+which leads, after a matrix inversion, spectral characteristics
+similar to schemes of higher order \citep{Lele_JCP1992}.
+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 \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 terms.
+advective fluxes \eqref{Eq_tra_adv} is kept at $2^{nd}$ order.
+The expression \textit{$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 the continuity equation,
+and the momentum advection and pressure 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 using
 \textit{cen4}. Furthermore, it must be used in conjunction with an explicit
 diffusion operator to produce a sensible solution. The time-stepping is also
 performed using a leapfrog scheme in conjunction with an Asselin time-filter,
 so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
 At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), an
 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 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})]
          {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}
+it is not non-diffusive, $i.e.$ the global variance of a tracer is not preserved using CEN4.
+Furthermore, it must be used in conjunction with an explicit diffusion operator
+to produce a sensible solution.
+As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction
+with an Asselin time-filter, so $T$ in (\ref{Eq_tra_adv_cen4}) is the \textit{now} tracer.
+At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface),
+an 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 a reduction to a second order scheme
+for these near boundary grid points.
+% -------------------------------------------------------------------------------------------------------------
+%        FCT scheme
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})]
+         {Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct}=true)}
 \label{TRA_adv_tvd}
+In the Total Variance Dissipation (TVD) formulation, the tracer at velocity
+points is evaluated using a combination of an upstream and a centred scheme.
+For example, in the $i$-direction :
+\begin{equation} \label{Eq_tra_adv_tvd}
+The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~=~\textit{true}.
+Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level)
+and vertical direction by setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
+FCT implementation can be found in the \mdl{traadv\_fct} module.
+In FCT formulation, the tracer at velocity 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_fct}
 \begin{split}
 \tau _u^{ups}&= \begin{cases}
 …
               \end{cases}     \\
 \\
 \tau _u^{tvd}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen2} -\tau _u^{ups} } \right)
+\tau _u^{fct}&=\tau _u^{ups} +c_u \;\left( {\tau _u^{cen} -\tau _u^{ups} } \right)
 \end{split}
 \end{equation}
 where $c_u$ is a flux limiter function taking values between 0 and 1.
+The FCT order is the one of the centred scheme used ($i.e.$ it depends on the setting of
+\np{nn\_fct\_h} and \np{nn\_fct\_v}.
 There exist many ways to define $c_u$, each corresponding to a different
+total variance decreasing scheme. The one chosen in \NEMO is described in
+\citet{Zalesak_JCP79}. $c_u$ only departs from $1$ when the advective term
+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{Levy_al_GRL01}; note that in this paper it is referred to as "FCT" (Flux corrected
+transport) rather than TVD. The TVD scheme is implemented in the \mdl{traadv\_tvd} module.
+For stability reasons (see \S\ref{STP}),
+$\tau _u^{cen2}$ is evaluated  in (\ref{Eq_tra_adv_tvd}) using the \textit{now} tracer while $\tau _u^{ups}$
+is evaluated using the \textit{before} tracer. In other words, the advective part of
+the scheme is time stepped with a leap-frog scheme while a forward scheme is
+used for the diffusive part.
+FCT scheme. The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
+$c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
+The resulting scheme is quite expensive but \emph{positive}.
+It can be used on both active and passive tracers.
+A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
+An additional option has been added controlled by \np{nn\_fct\_zts}. By setting this integer to
+a value larger than zero, a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction,
+but on the latter, a split-explicit time stepping is used, with a number of sub-timestep equals
+to \np{nn\_fct\_zts}. This option can be useful when the size of the timestep is limited
+by vertical advection \citep{Lemarie_OM2015}. Note that in this case, a similar split-explicit
+time stepping should be used on vertical advection of momentum to insure a better stability
+(see \S\ref{DYN_zad}).
+For stability reasons (see \S\ref{STP}), $\tau _u^{cen}$ is evaluated in (\ref{Eq_tra_adv_fct})
+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})]
+   {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{Levy_al_GRL01}. In its formulation, the tracer at velocity points
+\subsection[MUSCL scheme  (\np{ln\_traadv\_mus})]
+   {Monotone Upstream Scheme for Conservative Laws (MUSCL) (\np{ln\_traadv\_mus}=T)}
+\label{TRA_adv_mus}
+The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~=~\textit{true}.
+MUSCL implementation can be found in the \mdl{traadv\_mus} module.
+MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}. In its formulation, the tracer at velocity points
 is evaluated assuming a linear tracer variation between two $T$-points
 (Fig.\ref{Fig_adv_scheme}). For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_muscl}
+\begin{equation} \label{Eq_tra_adv_mus}
    \tau _u^{mus} = \left\{      \begin{aligned}
          &\tau _i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right)
 …
 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. The two MUSCL schemes are implemented
+in the \mdl{traadv\_tvd} and \mdl{traadv\_tvd2} modules.
+directed toward land, an upstream flux is used. This choice ensure
+the \textit{positive} character of the scheme.
+In addition, fluxes round a grid-point where a runoff is applied can optionally be
+computed using upstream fluxes (\np{ln\_mus\_ups}~=~\textit{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 the Cell
+Averaged QUICK scheme (Quadratic Upstream Interpolation for Convective
+Kinematics). For example, in the $i$-direction :
+The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~=~\textit{true}.
+UBS implementation can be found in the \mdl{traadv\_mus} module.
+The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme
+(Quadratic Upstream Interpolation for Convective Kinematics). It is an upstream-biased
+third order scheme based on an upstream-biased parabolic interpolation.
+For example, in the $i$-direction :
 \begin{equation} \label{Eq_tra_adv_ubs}
    \tau _u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{
 …
 This results in a dissipatively dominant (i.e. hyper-diffusive) truncation
 error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of the advection
 scheme is similar to that reported in \cite{Farrow1995}.
+error \citep{Shchepetkin_McWilliams_OM05}. The overall performance of
+ the advection scheme is similar to that reported in \cite{Farrow1995}.
 It is a relatively good compromise between accuracy and smoothness.
 It is not a \emph{positive} scheme, meaning that false extrema are permitted,
+Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted,
 but the amplitude of such are significantly reduced over the centred second
 order method. Nevertheless it is not recommended that it should be applied
 to a passive tracer that requires positivity.
+or fourth order method. therefore it is not recommended that it should be
+applied to a passive tracer that requires positivity.
 The intrinsic diffusion of UBS makes its use risky in the vertical direction
 where the control of artificial diapycnal fluxes is of paramount importance.
 Therefore the vertical flux is evaluated using the TVD scheme when
 \np{ln\_traadv\_ubs}=true.
+where the control of artificial diapycnal fluxes is of paramount importance \citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
+Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme
+or a $4^th$ order COMPACT scheme (\np{nn\_cen\_v}=2 or 4).
 For stability reasons  (see \S\ref{STP}),
 the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order centred scheme)
 is evaluated using the \textit{now} tracer (centred in time) while the
 second term (which is the diffusive part of the scheme), is
+the first term  in \eqref{Eq_tra_adv_ubs} (which corresponds to a second order
+centred scheme) is evaluated using the \textit{now} tracer (centred in time)
+while the second term (which is the diffusive part of the scheme), is
 evaluated using the \textit{before} tracer (forward in time).
 This choice is discussed by \citet{Webb_al_JAOT98} in the context of the
 …
 substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
+Four different options are possible for the vertical
+component used in the UBS scheme. $\tau _w^{ubs}$ can be evaluated
+using either \textit{(a)} a centred $2^{nd}$ order scheme, or  \textit{(b)}
+a TVD scheme, or  \textit{(c)} an interpolation based on conservative
+parabolic splines following the \citet{Shchepetkin_McWilliams_OM05}
+implementation of UBS in ROMS, or  \textit{(d)} a UBS. The $3^{rd}$ case
+has dispersion properties similar to an eighth-order accurate conventional scheme.
+The current reference version uses method b)
+Note that :
+(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. Such a technique has been
+implemented and validated in ORCA05 with 301 levels. It is not available
+in the current reference version.
+(2) It is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
+Note that it is straightforward to rewrite \eqref{Eq_tra_adv_ubs} as follows:
 \begin{equation} \label{Eq_traadv_ubs2}
 \tau _u^{ubs} = \tau _u^{cen4} + \frac{1}{12} \left\{
 …
 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 v3.4 still uses
+ \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}.
+ %%%
+ \gmcomment{the change in UBS scheme has to be done}
+ %%%
+ $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO uses
+the computationally more efficient formulation \eqref{Eq_tra_adv_ubs}.
 % -------------------------------------------------------------------------------------------------------------
 …
 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 the ULTIMATE QUICKEST
+is used when \np{ln\_traadv\_qck}~=~\textit{true}.
+QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
+QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST
 limiter \citep{Leonard1991}. It has been implemented in NEMO by G. Reffray
 (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
 …
 direction where the control of artificial diapycnal fluxes is of paramount importance.
 Therefore the vertical flux is evaluated using the CEN2 scheme.
+This no longer guarantees the positivity of the scheme. The use of TVD in the vertical
+direction (as for the UBS case) should be implemented to restore this property.
+% -------------------------------------------------------------------------------------------------------------
+%        PPM scheme
+% -------------------------------------------------------------------------------------------------------------
+\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)
+\sgacomment{reference?}
+is based on a quadradic piecewise construction. Like the QCK scheme, it is associated
+with the ULTIMATE QUICKEST limiter \citep{Leonard1991}. It has been implemented
+in \NEMO by G. Reffray (MERCATOR-ocean) but is not yet offered in the reference
+version 3.3.
+This no longer guarantees the positivity of the scheme.
+The use of FCT in the vertical direction (as for the UBS case) should be implemented
+to restore this property.
+%%%gmcomment   :  Cross term are missing in the current implementation....
 % ================================================================
 …
 %-------------------------------------------------------------------------------------------------------------
+Options are defined through the  \ngn{namtra\_ldf} namelist variables.
+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 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}.
+Options are defined through the \ngn{namtra\_ldf} namelist variables.
+They are regrouped in four items, allowing to specify
+$(i)$   the type of operator used (none, laplacian, bilaplacian),
+$(ii)$  the direction along which the operator acts (iso-level, horizontal, iso-neutral),
+$(iii)$ some specific options related to the rotated operators ($i.e.$ non-iso-level operator), and
+$(iv)$  the specification of eddy diffusivity coefficient (either constant or variable in space and time).
+Item $(iv)$ will be described in Chap.\ref{LDF} .
+The direction along which the operators act is defined through the slope between this direction and the iso-level surfaces.
+The slope is computed in the \mdl{ldfslp} module and will also be described in Chap.~\ref{LDF}.
 The lateral diffusion of tracers is evaluated using a forward scheme,
 $i.e.$ the tracers appearing in its expression are the \textit{before} tracers in time,
+except for the pure vertical component that appears when a rotation tensor
+is used. This latter term is solved implicitly together with the
+vertical diffusion term (see \S\ref{STP}).
+% -------------------------------------------------------------------------------------------------------------
+%        Iso-level laplacian operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap})]
+         {Iso-level laplacian operator (lap) (\np{ln\_traldf\_lap}=true) }
+\label{TRA_ldf_lap}
+A laplacian diffusion operator ($i.e.$ a harmonic operator) acting along the model
+surfaces is given by:
+except for the pure vertical component that appears when a rotation tensor is used.
+This latter component is solved implicitly together with the vertical diffusion term (see \S\ref{STP}).
+When \np{ln\_traldf\_msc}~=~\textit{true}, a Method of Stabilizing Correction is used in which
+the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}.
+% -------------------------------------------------------------------------------------------------------------
+%        Type of operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, \np{ln\_traldf\_blp})]
+              {Type of operator (\np{ln\_traldf\_NONE}, \np{ln\_traldf\_lap}, or \np{ln\_traldf\_blp} = true) }
+\label{TRA_ldf_op}
+Three operator options are proposed and, one and only one of them must be selected:
+\begin{description}
+\item [\np{ln\_traldf\_NONE}] = true : no operator selected, the lateral diffusive tendency will not be
+applied to the tracer equation. This option can be used when the selected advection scheme
+is diffusive enough (MUSCL scheme for example).
+\item [ \np{ln\_traldf\_lap}] = true : a laplacian operator is selected. This harmonic operator
+takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $,
+where the gradient operates along the selected direction (see \S\ref{TRA_ldf_dir}),
+and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see Chap.~\ref{LDF}).
+\item [\np{ln\_traldf\_blp}] = true : a bilaplacian operator is selected. This biharmonic operator
+takes the following expression:
+$\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$
+where the gradient operats along the selected direction,
+and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$  (see Chap.~\ref{LDF}).
+In the code, the bilaplacian operator is obtained by calling the laplacian twice.
+\end{description}
+Both laplacian and bilaplacian operators ensure the total tracer variance decrease.
+Their primary role is to provide strong dissipation at the smallest scale supported
+by the grid while minimizing the impact on the larger scale features.
+The main difference between the two operators is the scale selectiveness.
+The bilaplacian damping time ($i.e.$ its spin down time) scales like $\lambda^{-4}$
+for disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones),
+whereas the laplacian damping time scales only like $\lambda^{-2}$.
+% -------------------------------------------------------------------------------------------------------------
+%        Direction of action
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Direction of action (\np{ln\_traldf\_lev}, \np{ln\_traldf\_hor}, \np{ln\_traldf\_iso}, \np{ln\_traldf\_triad})]
+              {Direction of action (\np{ln\_traldf\_lev}, \textit{...\_hor}, \textit{...\_iso}, or \textit{...\_triad} = true) }
+\label{TRA_ldf_dir}
+The choice of a direction of action determines the form of operator used.
+The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane
+when iso-level option is used (\np{ln\_traldf\_lev}~=~\textit{true})
+or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{z}-coordinate
+(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \textit{true}).
+The associated code can be found in the \mdl{traldf\_lap\_blp} module.
+The operator is a rotated (re-entrant) laplacian when the direction along which it acts
+does not coincide with the iso-level surfaces,
+that is when standard or triad iso-neutral option is used (\np{ln\_traldf\_iso} or
+ \np{ln\_traldf\_triad} equals \textit{true}, see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.),
+or when a horizontal ($i.e.$ geopotential) operator is demanded in \textit{s}-coordinate
+(\np{ln\_traldf\_hor} and \np{ln\_sco} equal \textit{true})
+\footnote{In this case, the standard iso-neutral operator will be automatically selected}.
+In that case, a rotation is applied to the gradient(s) that appears in the operator
+so that diffusive fluxes acts on the three spatial direction.
+The resulting discret form of the three operators (one iso-level and two rotated one)
+is given in the next two sub-sections.
+% -------------------------------------------------------------------------------------------------------------
+%       iso-level operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso})]
+         {Iso-level (bi-)laplacian operator ( \np{ln\_traldf\_iso}) }
+\label{TRA_ldf_lev}
+The laplacian diffusion operator acting along the model (\textit{i,j})-surfaces is given by:
 \begin{equation} \label{Eq_tra_ldf_lap}
 D_T^{lT} =\frac{1}{b_tT} \left( \;
+D_t^{lT} =\frac{1}{b_t} \left( \;
    \delta _{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta _{i+1/2} [T] \right]
 + \delta _{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta _{j+1/2} [T] \right]  \;\right)
 \end{equation}
+where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells.
+It is implemented in the \mdl{traadv\_lap} module.
+This lateral operator is computed in \mdl{traldf\_lap}. It is a \emph{horizontal}
+operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate with
+or without partial steps, but is simply an iso-level operator in the $s$-coordinate.
+It is thus used when, in addition to \np{ln\_traldf\_lap}=true, we have
+\np{ln\_traldf\_level}=true or \np{ln\_traldf\_hor}=\np{ln\_zco}=true.
+where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells
+and where zero diffusive fluxes is assumed across solid boundaries,
+first (and third in bilaplacian case) horizontal tracer derivative are masked.
+It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module.
+The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap}
+in order to compute the iso-level bilaplacian operator.
+It is a \emph{horizontal} operator ($i.e.$ acting along geopotential surfaces) in the $z$-coordinate
+with or without partial steps, but is simply an iso-level operator in the $s$-coordinate.
+It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~=~\textit{true},
+we have \np{ln\_traldf\_lev}~=~\textit{true} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~=~\textit{true}.
 In both cases, it significantly contributes to diapycnal mixing.
 It is therefore not recommended.
+It is therefore never recommended, even when using it in the bilaplacian case.
 Note that in the partial step $z$-coordinate (\np{ln\_zps}=true), tracers in horizontally
 …
 described in \S\ref{TRA_zpshde}.
+% -------------------------------------------------------------------------------------------------------------
+%        Rotated laplacian operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Rotated laplacian operator (iso) (\np{ln\_traldf\_lap})]
+         {Rotated laplacian operator (iso) (\np{ln\_traldf\_lap}=true)}
+% -------------------------------------------------------------------------------------------------------------
+%         Rotated laplacian operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Standard and triad rotated (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad})]
+               {Standard and triad (bi-)laplacian operator (\mdl{traldf\_iso}, \mdl{traldf\_triad}))}
+\label{TRA_ldf_iso_triad}
+%&&    Standard rotated (bi-)laplacian operator
+%&& ----------------------------------------------
+\subsubsection   [Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})]
+                 {Standard rotated (bi-)laplacian operator (\mdl{traldf\_iso})}
 \label{TRA_ldf_iso}
+If the Griffies trad scheme is not employed
+(\np{ln\_traldf\_grif}=true; see App.\ref{sec:triad}) 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}
 \begin{split}
 …
 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{Guilyardi_al_CD01}. An alternative scheme
+developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
+background horizontal diffusion \citep{Guilyardi_al_CD01}.
+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}.
+%&&     Triad rotated (bi-)laplacian operator
+%&&  -------------------------------------------
+\subsubsection   [Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})]
+                 {Triad rotated (bi-)laplacian operator (\np{ln\_traldf\_triad})}
+\label{TRA_ldf_triad}
+If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}=true ; see App.\ref{sec:triad})
+An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
 is also available in \NEMO (\np{ln\_traldf\_grif}=true). A complete description of
 the algorithm is given in App.\ref{sec:triad}.
-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 (bilap) (\np{ln\_traldf\_bilap})]
-         {Iso-level bilaplacian operator (bilap) (\np{ln\_traldf\_bilap}=true)}
-\label{TRA_ldf_bilap}
 The lateral fourth order bilaplacian operator on tracers is obtained by
 applying (\ref{Eq_tra_ldf_lap}) twice. The operator requires an additional assumption
 on boundary conditions: both first and third derivative terms normal to the
+coast are set to zero. It is used when, in addition to \np{ln\_traldf\_bilap}=true,
+we have \np{ln\_traldf\_level}=true, or both \np{ln\_traldf\_hor}=true and
+\np{ln\_zco}=false. In both cases, it can contribute diapycnal mixing,
+although less than in the laplacian case. It is therefore not recommended.
+Note that in the code, the bilaplacian routine does not call the laplacian
+routine twice but is rather a separate routine that can be found in the
+\mdl{traldf\_bilap} module. This is due to the fact that we introduce the
+eddy diffusivity coefficient, A, in the operator as:
+$\nabla \cdot \nabla \left( {A\nabla \cdot \nabla T} \right)$,
+instead of
+$-\nabla \cdot a\nabla \left( {\nabla \cdot a\nabla T} \right)$
+where $a=\sqrt{|A|}$ and $A<0$. This was a mistake: both formulations
+ensure the total variance decrease, but the former requires a larger
+number of code-lines.
+% -------------------------------------------------------------------------------------------------------------
+%        Rotated bilaplacian operator
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap})]
+         {Rotated bilaplacian operator (bilapg) (\np{ln\_traldf\_bilap}=true)}
+\label{TRA_ldf_bilapg}
+coast are set to zero.
 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, normal to the bottom and normal to the surface are set to zero. It can be found in the
+\mdl{traldf\_bilapg}.
+It is used when, in addition to \np{ln\_traldf\_bilap}=true, we have
+\np{ln\_traldf\_iso}= .true, or both \np{ln\_traldf\_hor}=true and \np{ln\_zco}=true.
+This rotated bilaplacian operator has never been seriously
+tested. There are no guarantees that it is either free of bugs or correctly formulated.
+Moreover, the stability range of such an operator will be probably quite
+narrow, requiring a significantly smaller time-step than the one used with an
+unrotated operator.
+coast, normal to the bottom and normal to the surface are set to zero.
+%&&    Option for the rotated operators
+%&& ----------------------------------------------
+\subsubsection   [Option for the rotated operators]
+                 {Option for the rotated operators}
+\label{TRA_ldf_options}
+\np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators)
+\np{rn\_slpmax} = slope limit (both operators)
+\np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only)
+\np{rn\_sw\_triad} =1 switching triad ; =0 all 4 triads used (triad only)
+\np{ln\_botmix\_triad} = lateral mixing on bottom (triad only)
 % ================================================================
 …
 %--------------------------------------------------------------------------------------------------------------
 Options are defined through the  \ngn{namzdf} namelist variables.
+Options are defined through the \ngn{namzdf} namelist variables.
 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 thickness of the top model layer.
+Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components ($i.e.$ atmosphere, sea-ice, land),
+the change in the heat and salt content of the surface layer of the ocean is due both
+to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
+ and to the heat and salt content of the mass exchange.
+\sgacomment{ the following does not apply to the release to which this documentation is
+attached and so should not be included ....
+In a forthcoming release, these two parts, computed in the surface module (SBC), will be included directly
+in $Q_{ns}$, the surface heat flux and $F_{salt}$, the surface salt flux.
+The specification of these fluxes is further detailed in the SBC chapter (see \S\ref{SBC}).
+This change will provide a forcing formulation which is the same for any tracer (including temperature and salinity).
+In the current version, the situation is a little bit more complicated. }
+Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components
+($i.e.$ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer
+of the ocean is due both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$)
+and to the heat and salt content of the mass exchange. They are both included directly in $Q_{ns}$,
+the surface heat flux, and $F_{salt}$, the surface salt flux (see \S\ref{SBC} for further details).
+By doing this, the forcing formulation is the same for any tracer (including temperature and salinity).
 The surface module (\mdl{sbcmod}, see \S\ref{SBC}) provides the following
 …
 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface
 (i.e. the difference between the total surface heat flux and the fraction of the short wave flux that
+penetrates into the water column, see \S\ref{TRA_qsr})
+$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
+$\bullet$ $\textit{emp}_S$, an equivalent mass flux taking into account the effect of ice-ocean mass exchange
+$\bullet$ \textit{rnf}, the mass flux associated with runoff (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+The $\textit{emp}_S$ field is not simply the budget of evaporation-precipitation+freezing-melting because
+the sea-ice is not currently embedded in the ocean but levitates above it. There is no mass
+exchanged between the sea-ice and the ocean. Instead we only take into account the salt
+flux associated with the non-zero salinity of sea-ice, and the concentration/dilution effect
+due to the freezing/melting (F/M) process. These two parts of the forcing are then converted into
+an equivalent mass flux given by $\textit{emp}_S - \textit{emp}$. As a result of this mess,
+the surface boundary condition on temperature and salinity is applied as follows:
+In the nonlinear free surface case (\key{vvl} is defined):
+penetrates into the water column, see \S\ref{TRA_qsr}) plus the heat content associated with
+of the mass exchange with the atmosphere and lands.
+$\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...)
+$\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation)
+ and possibly with the sea-ice and ice-shelves.
+$\bullet$ \textit{rnf}, the mass flux associated with runoff
+(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
+$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt,
+(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
+The surface boundary condition on temperature and salinity is applied as follows:
 \begin{equation} \label{Eq_tra_sbc}
+\begin{aligned}
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^t  & \\
+& F^S =\frac{ 1 }{\rho _o  \,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\
+ \end{aligned}
+\end{equation}
+where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
+($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
+divergence of odd and even time step (see \S\ref{STP}).
+In the linear free surface case (\np{ln\_linssh}~=~\textit{true}),
+an additional term has to be added on both temperature and salinity.
+On temperature, this term remove the heat content associated with mass exchange
+that has been added to $Q_{ns}$. On salinity, this term mimics the concentration/dilution effect that
+would have resulted from a change in the volume of the first level.
+The resulting surface boundary condition is applied as follows:
+\begin{equation} \label{Eq_tra_sbc_lin}
 \begin{aligned}
  &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }
 …
+%
 & F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} }
            &\overline{ \left( (\textit{emp}_S - \textit{emp})\;\left. S \right|_{k=1}  \right) }^t   & \\
+           &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\
  \end{aligned}
 \end{equation}
+In the linear free surface case (\key{vvl} not defined):
+\begin{equation} \label{Eq_tra_sbc_lin}
+\begin{aligned}
+ &F^T = \frac{ 1 }{\rho _o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns} }^t  & \\
+%
+& F^S =\frac{ 1 }{\rho _o \,\left. e_{3t} \right|_{k=1} }
+           &\overline{ \left( \textit{emp}_S\;\left. S \right|_{k=1}  \right) }^t   & \\
+ \end{aligned}
+\end{equation}
+where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps
+($t-\rdt/2$ and $t+\rdt/2$). Such time averaging prevents the
+divergence of odd and even time step (see \S\ref{STP}).
+The two set of equations, \eqref{Eq_tra_sbc} and \eqref{Eq_tra_sbc_lin}, are obtained
+by assuming that the temperature of precipitation and evaporation are equal to
+the ocean surface temperature and that their salinity is zero. Therefore, the heat content
+of the \textit{emp} budget must be added to the temperature equation in the variable volume case,
+while it does not appear in the constant volume case. Similarly, the \textit{emp} budget affects
+the ocean surface salinity in the constant volume case (through the concentration dilution effect)
+while it does not appears explicitly in the variable volume case since salinity change will be
+induced by volume change. In both constant and variable volume cases, surface salinity
+will change with ice-ocean salt flux and F/M flux (both contained in $\textit{emp}_S - \textit{emp}$) without mass exchanges.
+Note that the concentration/dilution effect due to F/M is computed using
+a constant ice salinity as well as a constant ocean salinity.
+This approximation suppresses the correlation between \textit{SSS}
+and F/M flux, allowing the ice-ocean salt exchanges to be conservative.
+Indeed, if this approximation is not made, even if the F/M budget is zero
+on average over the whole ocean domain and over the seasonal cycle,
+the associated salt flux is not zero, since sea-surface salinity and F/M flux are
+intrinsically correlated (high \textit{SSS} are found where freezing is
+strong whilst low \textit{SSS} is usually associated with high melting areas).
+Even using this approximation, an exact conservation of heat and salt content
+is only achieved in the variable volume case. In the constant volume case,
+there is a small imbalance associated with the product $(\partial_t\eta - \textit{emp}) * \textit{SSS}$.
+Nevertheless, the salt content variation is quite small and will not induce
+a long term drift as there is no physical reason for $(\partial_t\eta - \textit{emp})$
+and \textit{SSS} to be correlated \citep{Roullet_Madec_JGR00}.
+Note that, while quite small, the imbalance in the constant volume case is larger
+Note that an exact conservation of heat and salt content is only achieved with non-linear free surface.
+In the linear free surface case, there is a small imbalance. The imbalance is larger
 than the imbalance associated with the Asselin time filter \citep{Leclair_Madec_OM09}.
 This is the reason why the modified filter is not applied in the constant volume case.
+This is the reason why the modified filter is not applied in the linear free surface case (see \S\ref{STP}).
 % -------------------------------------------------------------------------------------------------------------
 …
 ($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform
 chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb}
+in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation:
+(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed
+time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll
+by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB
+formulation is used to calculate both the phytoplankton light limitation in PISCES
+or LOBSTER and the oceanic heating rate.
+in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
+\begin{description}
+\item[\np{nn\_chdta}=0]
+a constant 0.05 g.Chl/L value everywhere ;
+\item[\np{nn\_chdta}=1]
+an observed time varying chlorophyll deduced from satellite surface ocean color measurement
+spread uniformly in the vertical direction ;
+\item[\np{nn\_chdta}=2]
+same as previous case except that a vertical profile of chlorophyl is used.
+Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
+\item[\np{ln\_qsr\_bio}=true]
+simulated time varying chlorophyll by TOP biogeochemical model.
+In this case, the RGB formulation is used to calculate both the phytoplankton
+light limitation in PISCES or LOBSTER and the oceanic heating rate.
+\end{description}
 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}.
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]     \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_Irradiance.pdf}
+\includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance}
 \caption{    \label{Fig_traqsr_irradiance}
 Penetration profile of the downward solar irradiance calculated by four models.
 …
 \label{TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
 \namdisplay{namtra_bbc}
+\namdisplay{nambbc}
 %--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]     \begin{center}
 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_TRA_geoth.pdf}
+\includegraphics[width=1.0\textwidth]{Fig_TRA_geoth}
 \caption{   \label{Fig_geothermal}
 Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!t]   \begin{center}
 \includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_BBL_adv.pdf}
+\includegraphics[width=0.7\textwidth]{Fig_BBL_adv}
 \caption{   \label{Fig_bbl}
 Advective/diffusive Bottom Boundary Layer. The BBL parameterisation is
 …
 \subsection[DMP\_TOOLS]{Generating resto.nc using DMP\_TOOLS}
+DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input. This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
+DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$.
+Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled
+and run on the same machine as the NEMO model. A mesh\_mask.nc file for the model configuration is required as an input.
+This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1.
+The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work.
+The \nl{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for the restoration coefficient.
 %--------------------------------------------nam_dmp_create-------------------------------------------------
 \namdisplay{nam_dmp_create}
+\namtools{namelist_dmp}
 %-------------------------------------------------------------------------------------------------------
 \np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and should be the same as specified in \nl{namcfg}. The variable \nl{lzoom} is used to specify that the damping is being used as in case \textit{a} above to provide boundary conditions to a zoom configuration. In the case of the arctic or antarctic zoom configurations this includes some specific treatment. Otherwise damping is applied to the 6 grid points along the ocean boundaries. The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in the \nl{nam\_zoom\_dmp} name list.
+The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations. \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for the ORCA4, ORCA2 and ORCA05 configurations. If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference configurations with previous model versions. \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. This option only has an effect if \np{ln\_full\_field} is true. \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. Finally \np{ln\_custom} specifies that the custom module will be called. This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
+The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}. Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to the full values of a 10$^{\circ}$ latitud band. This is often used because of the short adjustment time scale in the equatorial region \citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.
+The remaining switch namelist variables determine the spatial variation of the restoration coefficient in non-zoom configurations.
+\np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain.
+\np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea
+for the ORCA4, ORCA2 and ORCA05 configurations.
+If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as
+a function of the model number. This option is included to allow backwards compatability of the ORCA2 reference
+configurations with previous model versions.
+\np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines.
+This option only has an effect if \np{ln\_full\_field} is true.
+\np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer.
+Finally \np{ln\_custom} specifies that the custom module will be called.
+This module is contained in the file custom.F90 and can be edited by users. For example damping could be applied in a specific region.
+The restoration coefficient can be set to zero in equatorial regions by specifying a positive value of \np{nn\_hdmp}.
+Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to
+the full values of a 10\deg latitud band.
+This is often used because of the short adjustment time scale in the equatorial region
+\citep{Reverdin1991, Fujio1991, Marti_PhD92}. The time scale associated with the damping depends on the depth as a
+hyperbolic tangent, with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.
 % ================================================================
 …
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
 \subsection{Equation of State (\np{nn\_eos} = 0, 1 or 2)}
+\subsection{Equation Of Seawater (\np{nn\_eos} = -1, 0, or 1)}
 \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 seawater volumic mass, $\rho$,
+abusively called density, is a non linear empirical function of \textit{in situ}
+temperature, salinity and pressure. The reference equation of state is that
+defined by the Joint Panel on Oceanographic Tables and Standards
+\citep{UNESCO1983}. It was the standard equation of state used in early
+releases of OPA. However, even though this computation is fully vectorised,
+it is quite time consuming ($15$ to $20${\%} of the total CPU time) since
+it requires the prior computation of the \textit{in situ} temperature from the
+model \textit{potential} temperature using the \citep{Bryden1973} polynomial
+for adiabatic lapse rate and a $4^th$ order Runge-Kutta integration scheme.
+Since OPA6, we have used the \citet{JackMcD1995} equation of state for
+seawater instead. It allows the computation of the \textit{in situ} ocean density
+directly as a function of \textit{potential} temperature relative to the surface
+(an \NEMO variable), the practical salinity (another \NEMO variable) and the
+pressure (assuming no pressure variation along geopotential surfaces, $i.e.$
+the pressure in decibars is approximated by the depth in meters).
+Both the \citet{UNESCO1983} and \citet{JackMcD1995} equations of state
+have exactly the same except that the values of the various coefficients have
+been adjusted by \citet{JackMcD1995} in order to directly use the \textit{potential}
+temperature instead of the \textit{in situ} one. This reduces the CPU time of the
+\textit{in situ} density computation to about $3${\%} of the total CPU time,
+while maintaining a quite accurate equation of state.
+In the computer code, a \textit{true} density anomaly, $d_a= \rho / \rho_o - 1$,
+is computed, with $\rho_o$ a reference volumic mass. Called \textit{rau0}
+in the code, $\rho_o$ is defined in \mdl{phycst}, and a value of $1,035~Kg/m^3$.
+The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship
+linking seawater density, $\rho$, to a number of state variables,
+most typically temperature, salinity and pressure.
+Because density gradients control the pressure gradient force through the hydrostatic balance,
+the equation of state provides a fundamental bridge between the distribution of active tracers
+and the fluid dynamics. Nonlinearities of the EOS are of major importance, in particular
+influencing the circulation through determination of the static stability below the mixed layer,
+thus controlling rates of exchange between the atmosphere  and the ocean interior \citep{Roquet_JPO2015}.
+Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{UNESCO1983})
+or TEOS-10 \citep{TEOS10} standards should be used anytime a simulation of the real
+ocean circulation is attempted \citep{Roquet_JPO2015}.
+The use of TEOS-10 is highly recommended because
+\textit{(i)} it is the new official EOS,
+\textit{(ii)} it is more accurate, being based on an updated database of laboratory measurements, and
+\textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature
+and practical salinity for EOS-980, both variables being more suitable for use as model variables
+\citep{TEOS10, Graham_McDougall_JPO13}.
+EOS-80 is an obsolescent feature of the NEMO system, kept only for backward compatibility.
+For process studies, it is often convenient to use an approximation of the EOS. To that purposed,
+a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available.
+In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$,
+is computed, with $\rho_o$ a reference density. Called \textit{rau0}
+in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$.
 This is a sensible choice for the reference density used in a Boussinesq ocean
 climate model, as, with the exception of only a small percentage of the ocean,
+density in the World Ocean varies by no more than 2$\%$ from $1,035~kg/m^3$
+\citep{Gill1982}.
+Options are defined through the  \ngn{nameos} namelist variables.
+The default option (namelist parameter \np{nn\_eos}=0) is the \citet{JackMcD1995}
+equation of state. Its use is highly recommended. However, for process studies,
+it is often convenient to use a linear approximation of the density.
+density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}.
+Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos}
+which controls the EOS used (=-1 for TEOS10 ; =0 for EOS-80 ; =1 for S-EOS).
+\begin{description}
+\item[\np{nn\_eos}$=-1$] the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.
+The accuracy of this approximation is comparable to the TEOS-10 rational function approximation,
+but it is optimized for a boussinesq fluid and the polynomial expressions have simpler
+and more computationally efficient expressions for their derived quantities
+which make them more adapted for use in ocean models.
+Note that a slightly higher precision polynomial form is now used replacement of the TEOS-10
+rational function approximation for hydrographic data analysis  \citep{TEOS10}.
+A key point is that conservative state variables are used:
+Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$).
+The pressure in decibars is approximated by the depth in meters.
+With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. It is set to
+$C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
+Choosing polyTEOS10-bsq implies that the state variables used by the model are
+$\Theta$ and $S_A$. In particular, the initial state deined by the user have to be given as
+\textit{Conservative} Temperature and \textit{Absolute} Salinity.
+In addition, setting \np{ln\_useCT} to \textit{true} convert the Conservative SST to potential SST
+prior to either computing the air-sea and ice-sea fluxes (forced mode)
+or sending the SST field to the atmosphere (coupled mode).
+\item[\np{nn\_eos}$=0$] the polyEOS80-bsq equation of seawater is used.
+It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized
+to accurately fit EOS80 (Roquet, personal comm.). The state variables used in both the EOS80
+and the ocean model are:
+the Practical Salinity ((unit: psu, notation: $S_p$)) and Potential Temperature (unit: $^{\circ}C$, notation: $\theta$).
+The pressure in decibars is approximated by the depth in meters.
+With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature,
+salinity and pressure \citep{UNESCO1983}. Nevertheless, a severe assumption is made in order to
+have a heat content ($C_p T_p$) which is conserved by the model: $C_p$ is set to a constant
+value, the TEOS10 value.
+\item[\np{nn\_eos}$=1$] a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
+the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.)
+(see also \citet{Roquet_JPO2015}). It provides a simplistic linear representation of both
+cabbeling and thermobaricity effects which is enough for a proper treatment of the EOS
+in theoretical studies \citep{Roquet_JPO2015}.
 With such an equation of state there is no longer a distinction between
+\textit{in situ} and \textit{potential} density and both cabbeling and thermobaric
+effects are removed.
+Two linear formulations are available: a function of $T$ only (\np{nn\_eos}=1)
+and a function of both $T$ and $S$ (\np{nn\_eos}=2):
+\begin{equation} \label{Eq_tra_eos_linear}
+\textit{conservative} and \textit{potential} temperature, as well as between \textit{absolute}
+and \textit{practical} salinity.
+S-EOS takes the following expression:
+\begin{equation} \label{Eq_tra_S-EOS}
 \begin{split}
+  d_a(T)       &=  \rho (T)      /  \rho_o   - 1     =  \  0.0285         -  \alpha   \;T     \\
+  d_a(T,S)    &=  \rho (T,S)   /  \rho_o   - 1     =  \  \beta \; S       -  \alpha   \;T
+  d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\
+                                & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a  \\
+                                & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\
+  with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3
 \end{split}
 \end{equation}
+where $\alpha$ and $\beta$ are the thermal and haline expansion
+coefficients, and $\rho_o$, the reference volumic mass, $rau0$.
+($\alpha$ and $\beta$ can be modified through the \np{rn\_alpha} and
+\np{rn\_beta} namelist variables). Note that when $d_a$ is a function
+of $T$ only (\np{nn\_eos}=1), the salinity is a passive tracer and can be
+used as such.
+% -------------------------------------------------------------------------------------------------------------
+%        Brunt-Vais\"{a}l\"{a} Frequency
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
+where the computer name of the coefficients as well as their standard value are given in \ref{Tab_SEOS}.
+In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing
+the associated coefficients.
+Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS.
+setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS.
+Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S.
+\end{description}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{table}[!tb]
+\begin{center} \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|}
+\hline
+coeff.   & computer name   & S-EOS     &  description                      \\ \hline
+$a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline
+$b_0$       & \np{rn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline
+$\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline
+$\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline
+$\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline
+$\mu_1$     & \np{rn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline
+$\mu_2$     & \np{rn\_mu2}    & 1.1090 $10^{-5}$ &  thermobaric coeff. in S            \\ \hline
+\end{tabular}
+\caption{ \label{Tab_SEOS}
+Standard value of S-EOS coefficients. }
+\end{center}
+\end{table}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+% -------------------------------------------------------------------------------------------------------------
+%        Brunt-V\"{a}is\"{a}l\"{a} Frequency
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Brunt-V\"{a}is\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
 \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 for $N^2$
+ depends on the type of equation of state used (\np{nn\_eos} namelist parameter).
+For \np{nn\_eos}=0 (\citet{JackMcD1995} equation of state), the \citet{McDougall1987}
+polynomial expression is used (with the pressure in decibar approximated by
+the depth in meters):
+An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a}
+ frequency) is of paramount importance as determine the ocean stratification and
+ is used in several ocean parameterisations (namely TKE, GLS, Richardson number dependent
+ vertical diffusion, enhanced vertical diffusion, non-penetrative convection, tidal mixing
+ parameterisation, iso-neutral diffusion). In particular, $N^2$ has to be computed at the local pressure
+ (pressure in decibar being approximated by the depth in meters). The expression for $N^2$
+ is given by:
 \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$ and $\beta$ are the thermal and haline expansion coefficients.
-They are a function of  $\overline{T}^{\,k+1/2},\widetilde{S}=\overline{S}^{\,k+1/2} - 35.$,
-and  $z_w$, with $T$ the \textit{potential} temperature and $\widetilde{S}$ a salinity anomaly.
-Note that both $\alpha$ and $\beta$ depend on \textit{potential}
-temperature and salinity which are averaged at $w$-points prior
-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{nn\_eos}=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}.
+% -------------------------------------------------------------------------------------------------------------
+%        Specific Heat
+% -------------------------------------------------------------------------------------------------------------
+\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 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.
+where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS,
+and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients.
+The coefficients are a polynomial function of temperature, salinity and depth which expression
+depends on the chosen EOS. They are computed through \textit{eos\_rab}, a \textsc{Fortran}
+function that can be found in \mdl{eosbn2}.
 % -------------------------------------------------------------------------------------------------------------
 …
 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 freezing
 point is computed through \textit{tfreez}, a \textsc{Fortran} function that can be found
+point is computed through \textit{eos\_fzp}, a \textsc{Fortran} function that can be found
 in \mdl{eosbn2}.
+% -------------------------------------------------------------------------------------------------------------
+%        Potential Energy
+% -------------------------------------------------------------------------------------------------------------
+%\subsection{Potential Energy anomalies}
+%\label{TRA_bn2}
+%    =====>>>>> TO BE written
+%
 % ================================================================
 …
 \label{TRA_zpshde}
+\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}
+\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 cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=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 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
+The partial cell properties at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom.
+So, only the bottom interpolation is explained below.
 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
 …
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 \begin{figure}[!p]    \begin{center}
 \includegraphics[width=0.9\textwidth]{./TexFiles/Figures/Partial_step_scheme.pdf}
+\includegraphics[width=0.9\textwidth]{Partial_step_scheme}
 \caption{   \label{Fig_Partial_step_scheme}
 Discretisation of the horizontal difference and average of tracers in the $z$-partial
 …
 \gmcomment{gm :   this last remark has to be done}
 %%%
+\end{document}

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

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

Context Navigation

Changeset 7351 for branches/2016/dev_INGV_UKMO_2016/DOC/TexFiles/Chapters/Chap_TRA.tex

Legend:

branches/2016/dev_INGV_UKMO_2016/DOC/TexFiles/Chapters/Chap_TRA.tex

Download in other formats: