Changeset 6289 for trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

Timestamp:

2016-02-05T00:47:05+01:00 (8 years ago)

Author:

gm

Message:

#1673 DOC of the trunk - Update, see associated wiki page for description

File:

• trunk/DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (19 diffs)

trunk/DOC/TexFiles/Chapters/Chap_TRA.tex

@@ -36 +36 @@
 (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 
+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 
@@ -45 +45 @@
 
 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}).
 
@@ -54 +54 @@
 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{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
@@ -68 +68 @@
 %-------------------------------------------------------------------------------------------------------------
 
-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( 
@@ -171 +172 @@
 %        2nd and 4th order centred schemes
 % -------------------------------------------------------------------------------------------------------------
-\subsection [centred schemes (CEN) (\np{ln\_traadv\_cen})]
-            {centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
+\subsection [Centred schemes (CEN) (\np{ln\_traadv\_cen})]
+            {Centred schemes (CEN) (\np{ln\_traadv\_cen}=true)}
 \label{TRA_adv_cen}
 
@@ -278 +279 @@
 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 
+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}).
@@ -405 +406 @@
 Estimated Streaming Terms (QUICKEST) scheme proposed by \citet{Leonard1979} 
 is used when \np{ln\_traadv\_qck}~=~\textit{true}. 
-QUICKEST implementation can be found in the \mdl{traadv\_mus} module.
+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 
@@ -415 +416 @@
 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.
+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....
@@ -431 +433 @@
 %-------------------------------------------------------------------------------------------------------------
  
-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 
@@ -475 +543 @@
 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}
@@ -527 +599 @@
 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)
 
 % ================================================================
@@ -593 +653 @@
 %--------------------------------------------------------------------------------------------------------------
 
-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. 
@@ -651 +711 @@
 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 
@@ -669 +723 @@
 $\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)
+
+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} }   
@@ -692 +757 @@
 %
 & 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}).
 
 % -------------------------------------------------------------------------------------------------------------
@@ -849 +875 @@
 \label{TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
-\namdisplay{namtra_bbc}
+\namdisplay{nambbc}
 %--------------------------------------------------------------------------------------------------------------
 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -1093 +1119 @@
 \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$^{\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}.  
 
 % ================================================================
@@ -1271 +1318 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%        Brunt-Vais\"{a}l\"{a} Frequency
-% -------------------------------------------------------------------------------------------------------------
-\subsection{Brunt-Vais\"{a}l\"{a} Frequency (\np{nn\_eos} = 0, 1 or 2)}
+%        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}
+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 
@@ -1332 +1379 @@
 \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"}
+\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