Index: branches/UKMO/nemo_v3_6_STABLE_copy/DOC/TexFiles/Chapters/Chap_TRA.tex
===================================================================
--- branches/UKMO/nemo_v3_6_STABLE_copy/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 5780)
+++ branches/UKMO/nemo_v3_6_STABLE_copy/DOC/TexFiles/Chapters/Chap_TRA.tex	(revision 6436)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter 1 � Ocean Tracers (TRA)
+% Chapter 1 ——— Ocean Tracers (TRA)
 % ================================================================
 \chapter{Ocean Tracers (TRA)}
@@ -36,16 +36,14 @@
 (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}).
 
@@ -56,6 +54,6 @@
 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 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
@@ -125,7 +123,7 @@
 \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}, 
@@ -133,7 +131,8 @@
 
 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.
+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{nam\_traadv}} namelist, by 
@@ -146,8 +145,8 @@
 
 Note that 
-(1) cen2, cen4 and TVD schemes require an explicit diffusion 
+(1) cen2 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
+(2) cen2, MUSCL2, and UBS are not \textit{positive} schemes
 \footnote{negative values can appear in an initially strictly positive tracer field 
 which is advected}
@@ -189,5 +188,5 @@
 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 
+configuration only. This is an obsolescent feature as the recommended 
 advection scheme for the ORCA configuration is TVD (see  \S\ref{TRA_adv_tvd}).
 
@@ -196,41 +195,4 @@
 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:
-\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 
-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.  
-
-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.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -270,4 +232,12 @@
 used for the diffusive part. 
 
+An additional option has been added controlled by \np{ln\_traadv\_tvd\_zts}. 
+By setting this logical to true, a TVD scheme is used on both horizontal and vertical direction, 
+but on the latter, a split-explicit time stepping is used, with 5 sub-timesteps. 
+This option can be useful when the value 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 ensure a better stability (see \np{ln\_dynzad\_zts} in \S\ref{DYN_zad}).
+
+
 % -------------------------------------------------------------------------------------------------------------
 %        MUSCL scheme  
@@ -296,10 +266,12 @@
 
 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, 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.
+
+Note that when using np{ln\_traadv\_msc\_ups}~=~true in addition to \np{ln\_traadv\_muscl}=true, 
+the MUSCL fluxes are replaced by upstream fluxes in vicinity of river mouths.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -416,18 +388,4 @@
 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.
 
 % ================================================================
@@ -464,5 +422,5 @@
 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)
@@ -661,16 +619,10 @@
 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 
@@ -679,22 +631,37 @@
 $\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).\\
+
+In the non-linear free surface case (\key{vvl} is defined), 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 (\key{vvl} is \textit{not} defined), 
+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} }   
@@ -702,50 +669,11 @@
 %
 & 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}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -821,11 +749,19 @@
 ($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}. 
@@ -859,5 +795,5 @@
 \label{TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
-\namdisplay{namtra_bbc}
+\namdisplay{nambbc}
 %--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -1103,5 +1039,10 @@
 \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-------------------------------------------------
@@ -1111,7 +1052,23 @@
 \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$^{\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}.  
 
 % ================================================================
@@ -1167,115 +1124,138 @@
 %        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.
-
-% -------------------------------------------------------------------------------------------------------------
-%        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}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1298,5 +1278,5 @@
 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}.  
 
@@ -1308,11 +1288,13 @@
 \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 
+\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 for the hydrostatic pressure 
 gradient (\mdl{dynhpg} module) to be active. 
 \gmcomment{STEVEN from gm : question: not sure of  what -to be active- means}
+
 Before taking horizontal gradients between the tracers next to the bottom, a linear 
 interpolation in the vertical is used to approximate the deeper tracer as if it actually 
@@ -1390,2 +1372,14 @@
 \gmcomment{gm :   this last remark has to be done}
 %%%
+
+If under ice shelf seas opened (\np{ln\_isfcav}=true), the partial cell properties 
+at the top are computed in the same way as for the bottom. Some extra variables are, 
+however, computed to reduce the flow generated at the top and bottom if $z*$ coordinates activated.
+The extra variables calculated and used by \S\ref{DYN_hpg_isf} are:
+
+$\bullet$ $\overline{T}_k^{\,i+1/2}$ as described in \eqref{Eq_zps_hde}
+
+$\bullet$ $\delta _{i+1/2} Z_{T_k} = \widetilde {Z}^{\,i}_{T_k}-Z^{\,i}_{T_k}$ to compute 
+the pressure gradient correction term used by \eqref{Eq_dynhpg_sco} in \S\ref{DYN_hpg_isf},
+ with $\widetilde {Z}_{T_k}$ the depth of the point $\widetilde {T}_{k}$ in case of $z^*$ coordinates 
+(this term = 0 in z-coordinates)
