Changeset 6040 for branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex

Timestamp:

2015-12-14T09:23:38+01:00 (8 years ago)

Author:

gm

Message:

#1613: vvl by default : start to update the DOC for change in vvl, LDF and solvers

File:

• branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (23 diffs)

branches/2015/dev_r5836_NOC3_vvl_by_default/DOC/TexFiles/Chapters/Chap_TRA.tex

@@ -1 +1 @@
 % ================================================================
-% Chapter 1 � Ocean Tracers (TRA)
+% Chapter 1 ——— Ocean Tracers (TRA)
 % ================================================================
 \chapter{Ocean Tracers (TRA)}
@@ -40 +40 @@
 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}
-%%%
+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 
@@ -50 +48 @@
 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 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}.
+equation for output (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~=~true), as described in Chap.~\ref{MISC}.
 
 $\ $\newline    % force a new ligne
@@ -82 +80 @@
 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.
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -113 +111 @@
 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 
@@ -119 +117 @@
 $\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 
@@ -125 +123 @@
 \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 
+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}).
+
+The choice of an advection scheme 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. Details 
-of the advection schemes are given below. The choice of an advection scheme 
+\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 not be advected !
+
+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 
-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
+(1) CEN and FCT schemes require an explicit diffusion operator 
+while the other schemes are diffusive enough so that they do not necessarily require 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}
@@ -163 +164 @@
 
 % -------------------------------------------------------------------------------------------------------------
+%        2nd and 4th order centred schemes
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [$2^{nd}$ and $4^{th}$ order centred schemes (CEN) (\np{ln\_traadv\_cen})]
+         {$2^{nd}$ and $4^{th}$ order 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 
+
+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 :
@@ -176 +179 @@
 \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. 
+CEN2 is computed in the \mdl{traadv\_cen} module.
+
+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}
 
-% -------------------------------------------------------------------------------------------------------------
 %        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 velocity points as 
+a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 
+For example, in the $i$-direction:
 \begin{equation} \label{Eq_tra_adv_cen4}
 \tau _u^{cen4} 
@@ -211 +201 @@
 \end{equation}
 
-Strictly speaking, the cen4 scheme is not a $4^{th}$ order advection scheme 
+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 
+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.
@@ -235 +226 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%        TVD scheme  
-% -------------------------------------------------------------------------------------------------------------
-\subsection   [Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd})]
-         {Total Variance Dissipation scheme (TVD) (\np{ln\_traadv\_tvd}=true)}
+%        FCT scheme  
+% -------------------------------------------------------------------------------------------------------------
+\subsection   [$2^{nd}$ and $4^{th}$ Flux Corrected Transport schemes (FCT) (\np{ln\_traadv\_fct})]
+         {$2^{nd}$ and $4^{th}$ 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 
+In the Flux Corrected Transport 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}
+\begin{equation} \label{Eq_tra_adv_fct}
 \begin{split}
 \tau _u^{ups}&= \begin{cases}
@@ -251 +242 @@
               \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}
@@ -260 +251 @@
 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.
+This scheme is tested and compared with MUSCL and a MPDATA scheme in \citet{Levy_al_GRL01}. 
+The FCT scheme is implemented in the \mdl{traadv\_fct} 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. 
+$\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}
+\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) has been 
@@ -281 +271 @@
 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)
@@ -296 +286 @@
 
 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. 
 
 % -------------------------------------------------------------------------------------------------------------
@@ -310 +296 @@
 \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 
+The UBS advection scheme (also often called UP3) 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 :
 \begin{equation} \label{Eq_tra_adv_ubs}
@@ -324 +310 @@
 
 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, 
 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. 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. 
-Therefore the vertical flux is evaluated using the TVD scheme when 
-\np{ln\_traadv\_ubs}=true.
+Therefore the vertical flux is evaluated using either a 2nd order FCT scheme 
+or a 4th 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 
@@ -350 +336 @@
 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 
+a FCT 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)
+The current reference version uses method (b).
+
+???
 
 Note that :
@@ -390 +380 @@
 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 
+ $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. Note the current version of NEMO still uses 
  \eqref{Eq_tra_adv_ubs}, not \eqref{Eq_traadv_ubs2}.
  %%%
@@ -416 +406 @@
 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.
 
 % ================================================================
@@ -1167 +1143 @@
 %        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: $\degres C$, 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}\,\degres 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: $\degres 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.
+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{nn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline
+$b_0$       & \np{nn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline
+$\lambda_1$ & \np{nn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline
+$\lambda_2$ & \np{nn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline
+$\nu$       & \np{nn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline
+$\mu_1$     & \np{nn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline
+$\mu_2$     & \np{nn\_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}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1232 +1263 @@
 
 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): 
+ 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}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1298 +1297 @@
 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
+%
+
 
 % ================================================================