Changeset 11422 for NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_TRA.tex

Timestamp:

2019-08-08T15:40:47+02:00 (5 years ago)

Author:

jchanut

Message:

#1791, merge with trunk

Location:

NEMO/branches/2019/fix_vvl_ticket1791/doc

Files:

• . (modified) (1 prop)
• latex (modified) (1 prop)
• latex/NEMO (modified) (1 prop)
• latex/NEMO/subfiles/chap_TRA.tex (modified) (70 diffs)

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures

NEMO/branches/2019/fix_vvl_ticket1791/doc/latex/NEMO/subfiles/chap_TRA.tex

@@ -55 +55 @@
 
 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}~\forcode{= .true.}), as described in \autoref{chap:DIA}.
+(\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}.
 
 % ================================================================
 % Tracer Advection
 % ================================================================
-\section{Tracer advection (\protect\mdl{traadv})}
+\section[Tracer advection (\textit{traadv.F90})]
+{Tracer advection (\protect\mdl{traadv})}
 \label{sec:TRA_adv}
 %------------------------------------------namtra_adv-----------------------------------------------------
@@ -81 +82 @@
 Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ 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, \ie \np{ln\_linssh}~\forcode{= .true.}).
+(which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie \np{ln\_linssh}\forcode{ = .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
@@ -90 +91 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[]{Fig_adv_scheme}
+    \includegraphics[width=\textwidth]{Fig_adv_scheme}
     \caption{
       \protect\label{fig:adv_scheme}
@@ -119 +120 @@
 \begin{description}
 \item[linear free surface:]
-  (\np{ln\_linssh}~\forcode{= .true.})
+  (\np{ln\_linssh}\forcode{ = .true.})
   the first level thickness is constant in time:
   the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on
@@ -127 +128 @@
   the first level tracer value.
 \item[non-linear free surface:]
-  (\np{ln\_linssh}~\forcode{= .false.})
+  (\np{ln\_linssh}\forcode{ = .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 fluxes through the surface are also zero.
@@ -136 +137 @@
 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 \citep{Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}.
-
-The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) is
+two quantities that are not correlated \citep{roullet.madec_JGR00, griffies.pacanowski.ea_MWR01, campin.adcroft.ea_OM04}.
+
+The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco?}) is
 the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity
 (see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or
@@ -183 +184 @@
 %        2nd and 4th order centred schemes
 % -------------------------------------------------------------------------------------------------------------
-\subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}~\forcode{= .true.})}
+\subsection[CEN: Centred scheme (\forcode{ln_traadv_cen = .true.})]
+{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}\forcode{ = .true.})}
 \label{subsec:TRA_adv_cen}
 
 %        2nd order centred scheme  
 
-The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~\forcode{= .true.}.
+The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{ = .true.}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$.
@@ -220 +222 @@
   \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2}
 \end{equation}
-In the vertical direction (\np{nn\_cen\_v}~\forcode{= 4}),
-a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}.
+In the vertical direction (\np{nn\_cen\_v}\forcode{ = 4}),
+a $4^{th}$ COMPACT interpolation has been prefered \citep{demange_phd14}.
 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion,
-spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 
+spectral characteristics similar to schemes of higher order \citep{lele_JCP92}. 
 
 Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but
@@ -250 +252 @@
 %        FCT scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}~\forcode{= .true.})}
+\subsection[FCT: Flux Corrected Transport scheme (\forcode{ln_traadv_fct = .true.})]
+{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}\forcode{ = .true.})}
 \label{subsec:TRA_adv_tvd}
 
-The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~\forcode{= .true.}.
+The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}.
 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by
 setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$.
@@ -277 +280 @@
 (\ie it depends on the setting of \np{nn\_fct\_h} and \np{nn\_fct\_v}).
 There exist many ways to define $c_u$, each corresponding to a different FCT scheme.
-The one chosen in \NEMO is described in \citet{Zalesak_JCP79}.
+The one chosen in \NEMO is described in \citet{zalesak_JCP79}.
 $c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field.
 The resulting scheme is quite expensive but \textit{positive}.
 It can be used on both active and passive tracers.
-A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}.
+A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{levy.estublier.ea_GRL01}.
 
 An additional option has been added controlled by \np{nn\_fct\_zts}.
@@ -287 +290 @@
 a $2^{nd}$ order FCT scheme is used on both horizontal and vertical direction, but on the latter,
 a split-explicit time stepping is used, with a number of sub-timestep equals to \np{nn\_fct\_zts}.
-This option can be useful when the size of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.
+This option can be useful when the size of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}.
 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 \autoref{subsec:DYN_zad}).
@@ -300 +303 @@
 %        MUSCL scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}~\forcode{= .true.})}
+\subsection[MUSCL: Monotone Upstream Scheme for Conservative Laws (\forcode{ln_traadv_mus = .true.})]
+{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}\forcode{ = .true.})}
 \label{subsec:TRA_adv_mus}
 
-The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~\forcode{= .true.}.
+The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}\forcode{ = .true.}.
 MUSCL implementation can be found in the \mdl{traadv\_mus} module.
 
-MUSCL has been first implemented in \NEMO by \citet{Levy_al_GRL01}.
+MUSCL has been first implemented in \NEMO by \citet{levy.estublier.ea_GRL01}.
 In its formulation, the tracer at velocity points is evaluated assuming a linear tracer variation between
 two $T$-points (\autoref{fig:adv_scheme}).
@@ -331 +335 @@
 This choice ensure the \textit{positive} character of the scheme.
 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes
-(\np{ln\_mus\_ups}~\forcode{= .true.}).
+(\np{ln\_mus\_ups}\forcode{ = .true.}).
 
 % -------------------------------------------------------------------------------------------------------------
 %        UBS scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}~\forcode{= .true.})}
+\subsection[UBS a.k.a. UP3: Upstream-Biased Scheme (\forcode{ln_traadv_ubs = .true.})]
+{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})}
 \label{subsec:TRA_adv_ubs}
 
-The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~\forcode{= .true.}.
+The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}.
 UBS implementation can be found in the \mdl{traadv\_mus} module.
 
@@ -358 +363 @@
 
 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}.
+\citep{shchepetkin.mcwilliams_OM05}.
+The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}.
 It is a relatively good compromise between accuracy and smoothness.
 Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted,
@@ -367 +372 @@
 The intrinsic diffusion of UBS makes its use risky in the vertical direction where
 the control of artificial diapycnal fluxes is of paramount importance
-\citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}.
+\citep{shchepetkin.mcwilliams_OM05, demange_phd14}.
 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme
-(\np{nn\_cen\_v}~\forcode{= 2 or 4}).
+(\np{nn\_cen\_v}\forcode{ = 2 or 4}).
 
 For stability reasons (see \autoref{chap:STP}), the first term  in \autoref{eq:tra_adv_ubs}
@@ -376 +381 @@
 (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 QUICK advection scheme.
+This choice is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the QUICK advection scheme.
 UBS and QUICK schemes only differ by one coefficient.
-Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} leads to the QUICK advection scheme \citep{Webb_al_JAOT98}.
+Replacing 1/6 with 1/8 in \autoref{eq:tra_adv_ubs} leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}.
 This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded.
 Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme.
@@ -408 +413 @@
 %        QCK scheme  
 % -------------------------------------------------------------------------------------------------------------
-\subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}~\forcode{= .true.})}
+\subsection[QCK: QuiCKest scheme (\forcode{ln_traadv_qck = .true.})]
+{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}\forcode{ = .true.})}
 \label{subsec:TRA_adv_qck}
 
 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme
-proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}.
+proposed by \citet{leonard_CMAME79} is used when \np{ln\_traadv\_qck}\forcode{ = .true.}.
 QUICKEST implementation can be found in the \mdl{traadv\_qck} module.
 
 QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST limiter
-\citep{Leonard1991}.
+\citep{leonard_CMAME91}.
 It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module.
 The resulting scheme is quite expensive but \textit{positive}.
@@ -431 +437 @@
 % Tracer Lateral Diffusion
 % ================================================================
-\section{Tracer lateral diffusion (\protect\mdl{traldf})}
+\section[Tracer lateral diffusion (\textit{traldf.F90})]
+{Tracer lateral diffusion (\protect\mdl{traldf})}
 \label{sec:TRA_ldf}
 %-----------------------------------------nam_traldf------------------------------------------------------
@@ -453 +460 @@
 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 \autoref{chap:STP}).
-When \np{ln\_traldf\_msc}~\forcode{= .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}.
+When \np{ln\_traldf\_msc}\forcode{ = .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.debreu.ea_OM12}.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Type of operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})]{Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 
+\subsection[Type of operator (\texttt{ln\_traldf}\{\texttt{\_NONE,\_lap,\_blp}\})]
+{Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) } 
 \label{subsec:TRA_ldf_op}
 
@@ -465 +473 @@
 
 \begin{description}
-\item[\np{ln\_traldf\_NONE}~\forcode{= .true.}:]
+\item[\np{ln\_traldf\_NONE}\forcode{ = .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}~\forcode{= .true.}:]
+\item[\np{ln\_traldf\_lap}\forcode{ = .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 \autoref{subsec:TRA_ldf_dir}),
   and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}).
-\item[\np{ln\_traldf\_blp}~\forcode{= .true.}]:
+\item[\np{ln\_traldf\_blp}\forcode{ = .true.}]:
   a bilaplacian operator is selected.
   This biharmonic operator takes the following expression:
@@ -493 +501 @@
 %        Direction of action
 % -------------------------------------------------------------------------------------------------------------
-\subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})]{Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 
+\subsection[Action direction (\texttt{ln\_traldf}\{\texttt{\_lev,\_hor,\_iso,\_triad}\})]
+{Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) } 
 \label{subsec: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}~\forcode{= .true.}) or
+iso-level option is used (\np{ln\_traldf\_lev}\forcode{ = .true.}) or
 when a horizontal (\ie geopotential) operator is demanded in \textit{z}-coordinate
 (\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}).
@@ -519 +528 @@
 %       iso-level operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) }
+\subsection[Iso-level (bi-)laplacian operator (\texttt{ln\_traldf\_iso})]
+{Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso})}
 \label{subsec:TRA_ldf_lev}
 
@@ -537 +547 @@
 It is a \textit{horizontal} operator (\ie 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}~\forcode{= .true.},
-we have \np{ln\_traldf\_lev}~\forcode{= .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~\forcode{= .true.}.
+It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{ = .true.},
+we have \np{ln\_traldf\_lev}\forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{ = .true.}.
 In both cases, it significantly contributes to diapycnal mixing.
 It is therefore never recommended, even when using it in the bilaplacian case.
 
-Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}),
+Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}),
 tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom.
 In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level require a specific treatment.
@@ -550 +560 @@
 %         Rotated laplacian operator
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Standard and triad (bi -)laplacian operator}
+\subsection{Standard and triad (bi-)laplacian operator}
 \label{subsec:TRA_ldf_iso_triad}
 
-%&&    Standard rotated (bi -)laplacian operator
+%&&    Standard rotated (bi-)laplacian operator
 %&& ----------------------------------------------
-\subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})}
+\subsubsection[Standard rotated (bi-)laplacian operator (\textit{traldf\_iso.F90})]
+{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})}
 \label{subsec:TRA_ldf_iso}
 The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:PE_zdf})
@@ -574 +585 @@
 $r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and
 the surface along which the diffusion operator acts (\ie horizontal or iso-neutral surfaces).
-It is thus used when, in addition to \np{ln\_traldf\_lap}~\forcode{= .true.},
-we have \np{ln\_traldf\_iso}~\forcode{= .true.},
-or both \np{ln\_traldf\_hor}~\forcode{= .true.} and \np{ln\_zco}~\forcode{= .true.}.
+It is thus used when, in addition to \np{ln\_traldf\_lap}\forcode{ = .true.},
+we have \np{ln\_traldf\_iso}\forcode{ = .true.},
+or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}.
 The way these slopes are evaluated is given in \autoref{sec:LDF_slp}.
 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using
@@ -590 +601 @@
 This formulation conserves the tracer but does not ensure the decrease of the tracer variance.
 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without
-any additional background horizontal diffusion \citep{Guilyardi_al_CD01}.
-
-Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}),
+any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}.
+
+Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}),
 the horizontal derivatives at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment.
 They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}.
 
-%&&     Triad rotated (bi -)laplacian operator
+%&&     Triad rotated (bi-)laplacian operator
 %&&  -------------------------------------------
-\subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})}
+\subsubsection[Triad rotated (bi-)laplacian operator (\textit{ln\_traldf\_triad})]
+{Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})}
 \label{subsec:TRA_ldf_triad}
 
-If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad})
-
-An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases
-is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}).
+If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:triad})
+
+An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which ensures tracer variance decreases
+is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}).
 A complete description of the algorithm is given in \autoref{apdx:triad}.
 
@@ -632 +644 @@
 % Tracer Vertical Diffusion
 % ================================================================
-\section{Tracer vertical diffusion (\protect\mdl{trazdf})}
+\section[Tracer vertical diffusion (\textit{trazdf.F90})]
+{Tracer vertical diffusion (\protect\mdl{trazdf})}
 \label{sec:TRA_zdf}
 %--------------------------------------------namzdf---------------------------------------------------------
@@ -663 +676 @@
 
 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that
-in the case of explicit time stepping (\np{ln\_zdfexp}~\forcode{= .true.})
+in the case of explicit time stepping (\np{ln\_zdfexp}\forcode{ = .true.})
 there would be too restrictive a constraint on the time step.
 Therefore, the default implicit time stepping is preferred for the vertical diffusion since
 it overcomes the stability constraint.
-A forward time differencing scheme (\np{ln\_zdfexp}~\forcode{= .true.}) using
+A forward time differencing scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using
 a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative.
 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.
@@ -680 +693 @@
 %        surface boundary condition
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Surface boundary condition (\protect\mdl{trasbc})}
+\subsection[Surface boundary condition (\textit{trasbc.F90})]
+{Surface boundary condition (\protect\mdl{trasbc})}
 \label{subsec:TRA_sbc}
 
@@ -730 +744 @@
 Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:STP}).
 
-In the linear free surface case (\np{ln\_linssh}~\forcode{= .true.}), an additional term has to be added on
+In the linear free surface case (\np{ln\_linssh}\forcode{ = .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}$.
@@ -747 +761 @@
 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}.
+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 linear free surface case (see \autoref{chap:STP}).
 
@@ -753 +767 @@
 %        Solar Radiation Penetration 
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Solar radiation penetration (\protect\mdl{traqsr})}
+\subsection[Solar radiation penetration (\textit{traqsr.F90})]
+{Solar radiation penetration (\protect\mdl{traqsr})}
 \label{subsec:TRA_qsr}
 %--------------------------------------------namqsr--------------------------------------------------------
@@ -761 +776 @@
 
 Options are defined through the \ngn{namtra\_qsr} namelist variables.
-When the penetrative solar radiation option is used (\np{ln\_flxqsr}~\forcode{= .true.}),
+When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}),
 the solar radiation penetrates the top few tens of meters of the ocean.
-If it is not used (\np{ln\_flxqsr}~\forcode{= .false.}) all the heat flux is absorbed in the first ocean level.
+If it is not used (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level.
 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:PE_tra_T} and
 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 
@@ -792 +807 @@
 larger depths where it contributes to local heating.
 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen.
-In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}~\forcode{= .true.})
+In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.})
 a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths,
-leading to the following expression \citep{Paulson1977}:
+leading to the following expression \citep{paulson.simpson_JPO77}:
 \[
   % \label{eq:traqsr_iradiance}
@@ -805 +820 @@
 
 Such assumptions have been shown to provide a very crude and simplistic representation of
-observed light penetration profiles (\cite{Morel_JGR88}, see also \autoref{fig:traqsr_irradiance}).
+observed light penetration profiles (\cite{morel_JGR88}, see also \autoref{fig:traqsr_irradiance}).
 Light absorption in the ocean depends on particle concentration and is spectrally selective.
-\cite{Morel_JGR88} has shown that an accurate representation of light penetration can be provided by
+\cite{morel_JGR88} has shown that an accurate representation of light penetration can be provided by
 a 61 waveband formulation.
 Unfortunately, such a model is very computationally expensive.
-Thus, \cite{Lengaigne_al_CD07} have constructed a simplified version of this formulation in which
+Thus, \cite{lengaigne.menkes.ea_CD07} have constructed a simplified version of this formulation in which
 visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700nm).
 For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to the coefficients computed from
-the full spectral model of \cite{Morel_JGR88} (as modified by \cite{Morel_Maritorena_JGR01}),
+the full spectral model of \cite{morel_JGR88} (as modified by \cite{morel.maritorena_JGR01}),
 assuming the same power-law relationship.
 As shown in \autoref{fig:traqsr_irradiance}, this formulation, called RGB (Red-Green-Blue),
@@ -820 +835 @@
 The 2-bands formulation does not reproduce the full model very well.
 
-The RGB formulation is used when \np{ln\_qsr\_rgb}~\forcode{= .true.}.
+The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{ = .true.}.
 The RGB attenuation coefficients (\ie the inverses of the extinction length scales) are tabulated over
 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L
@@ -827 +842 @@
 
 \begin{description}
-\item[\np{nn\_chdta}~\forcode{= 0}]
+\item[\np{nn\_chdta}\forcode{ = 0}]
   a constant 0.05 g.Chl/L value everywhere ; 
-\item[\np{nn\_chdta}~\forcode{= 1}]
+\item[\np{nn\_chdta}\forcode{ = 1}]
   an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in
   the vertical direction;
-\item[\np{nn\_chdta}~\forcode{= 2}]
+\item[\np{nn\_chdta}\forcode{ = 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}~\forcode{= .true.}]
+  Following \cite{morel.berthon_LO89}, the profile is computed from the local surface chlorophyll value;
+\item[\np{ln\_qsr\_bio}\forcode{ = .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
@@ -856 +871 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[]{Fig_TRA_Irradiance}
+    \includegraphics[width=\textwidth]{Fig_TRA_Irradiance}
     \caption{
       \protect\label{fig:traqsr_irradiance}
@@ -865 +880 @@
       61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of
       (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$.
-      From \citet{Lengaigne_al_CD07}.
+      From \citet{lengaigne.menkes.ea_CD07}.
     }
   \end{center}
@@ -874 +889 @@
 %        Bottom Boundary Condition
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Bottom boundary condition (\protect\mdl{trabbc})}
+\subsection[Bottom boundary condition (\textit{trabbc.F90})]
+{Bottom boundary condition (\protect\mdl{trabbc})}
 \label{subsec:TRA_bbc}
 %--------------------------------------------nambbc--------------------------------------------------------
@@ -883 +899 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[]{Fig_TRA_geoth}
+    \includegraphics[width=\textwidth]{Fig_TRA_geoth}
     \caption{
       \protect\label{fig:geothermal}
-      Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{Emile-Geay_Madec_OS09}.
-      It is inferred from the age of the sea floor and the formulae of \citet{Stein_Stein_Nat92}.
+      Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{emile-geay.madec_OS09}.
+      It is inferred from the age of the sea floor and the formulae of \citet{stein.stein_N92}.
     }
   \end{center}
@@ -897 +913 @@
 This is the default option in \NEMO, and it is implemented using the masking technique.
 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling.
-This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{Stein_Stein_Nat92}),
+This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{stein.stein_N92}),
 but it warms systematically the ocean and acts on the densest water masses.
 Taking this flux into account in a global ocean model increases the deepest overturning cell
-(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}.
+(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{emile-geay.madec_OS09}.
 
 Options are defined through the  \ngn{namtra\_bbc} namelist variables.
@@ -907 +923 @@
 the \np{nn\_geoflx\_cst}, which is also a namelist parameter.
 When \np{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is introduced which is provided in
-the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:geothermal}) \citep{Emile-Geay_Madec_OS09}.
+the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:geothermal}) \citep{emile-geay.madec_OS09}.
 
 % ================================================================
 % Bottom Boundary Layer
 % ================================================================
-\section{Bottom boundary layer (\protect\mdl{trabbl} - \protect\key{trabbl})}
+\section[Bottom boundary layer (\textit{trabbl.F90} - \texttt{\textbf{key\_trabbl}})]
+{Bottom boundary layer (\protect\mdl{trabbl} - \protect\key{trabbl})}
 \label{sec:TRA_bbl}
 %--------------------------------------------nambbl---------------------------------------------------------
@@ -931 +948 @@
 sometimes over a thickness much larger than the thickness of the observed gravity plume.
 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of
-a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved.
-
-The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997},
+a sill \citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved.
+
+The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{beckmann.doscher_JPO97},
 is to allow a direct communication between two adjacent bottom cells at different levels,
 whenever the densest water is located above the less dense water.
@@ -939 +956 @@
 In the current implementation of the BBL, only the tracers are modified, not the velocities.
 Furthermore, it only connects ocean bottom cells, and therefore does not include all the improvements introduced by
-\citet{Campin_Goosse_Tel99}.
+\citet{campin.goosse_T99}.
 
 % -------------------------------------------------------------------------------------------------------------
 %        Diffusive BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}~\forcode{= 1})}
+\subsection[Diffusive bottom boundary layer (\forcode{nn_bbl_ldf = 1})]
+{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}\forcode{ = 1})}
 \label{subsec:TRA_bbl_diff}
 
@@ -955 +973 @@
 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and
 $A_l^\sigma$ the lateral diffusivity in the BBL.
-Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence,
+Following \citet{beckmann.doscher_JPO97}, the latter is prescribed with a spatial dependence,
 \ie in the conditional form
 \begin{equation}
@@ -983 +1001 @@
 %        Advective BBL
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}~\forcode{= 1..2})}
+\subsection[Advective bottom boundary layer (\forcode{nn_bbl_adv = [12]})]
+{Advective bottom boundary layer (\protect\np{nn\_bbl\_adv}\forcode{ = [12]})}
 \label{subsec:TRA_bbl_adv}
 
@@ -994 +1013 @@
 \begin{figure}[!t]
   \begin{center}
-    \includegraphics[]{Fig_BBL_adv}
+    \includegraphics[width=\textwidth]{Fig_BBL_adv}
     \caption{
       \protect\label{fig:bbl}
@@ -1014 +1033 @@
 %%%gmcomment   :  this section has to be really written
 
-When applying an advective BBL (\np{nn\_bbl\_adv}~\forcode{= 1..2}), an overturning circulation is added which
+When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{ = 1..2}), an overturning circulation is added which
 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope.
 The density difference causes dense water to move down the slope.
 
-\np{nn\_bbl\_adv}~\forcode{= 1}:
+\np{nn\_bbl\_adv}\forcode{ = 1}:
 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step
-(see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}.
+(see black arrow in \autoref{fig:bbl}) \citep{beckmann.doscher_JPO97}.
 It is a \textit{conditional advection}, that is, advection is allowed only
 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and
 if the velocity is directed towards greater depth (\ie $\vect U \cdot \nabla H > 0$).
 
-\np{nn\_bbl\_adv}~\forcode{= 2}:
+\np{nn\_bbl\_adv}\forcode{ = 2}:
 the downslope velocity is chosen to be proportional to $\Delta \rho$,
-the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}.
+the density difference between the higher cell and lower cell densities \citep{campin.goosse_T99}.
 The advection is allowed only  if dense water overlies less dense water on the slope
 (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$).
@@ -1041 +1060 @@
 The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity,
 and because no direct estimation of this parameter is available, a uniform value has been assumed.
-The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}.
+The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}.
 
 Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme.
@@ -1074 +1093 @@
 % Tracer damping
 % ================================================================
-\section{Tracer damping (\protect\mdl{tradmp})}
+\section[Tracer damping (\textit{tradmp.F90})]
+{Tracer damping (\protect\mdl{tradmp})}
 \label{sec:TRA_dmp}
 %--------------------------------------------namtra_dmp-------------------------------------------------
@@ -1109 +1129 @@
 In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) whereas
 it is zero in the interior of the model domain.
-The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}.
+The second case corresponds to the use of the robust diagnostic method \citep{sarmiento.bryan_JGR82}.
 It allows us to find the velocity field consistent with the model dynamics whilst
 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$).
@@ -1121 +1141 @@
 only below the mixed layer (defined either on a density or $S_o$ criterion).
 It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here
-\citep{Madec_al_JPO96}.
+\citep{madec.delecluse.ea_JPO96}.
 
 For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under
@@ -1129 +1149 @@
 % Tracer time evolution
 % ================================================================
-\section{Tracer time evolution (\protect\mdl{tranxt})}
+\section[Tracer time evolution (\textit{tranxt.F90})]
+{Tracer time evolution (\protect\mdl{tranxt})}
 \label{sec:TRA_nxt}
 %--------------------------------------------namdom-----------------------------------------------------
@@ -1137 +1158 @@
 
 Options are defined through the \ngn{namdom} namelist variables.
-The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09},
+The general framework for tracer time stepping is a modified leap-frog scheme \citep{leclair.madec_OM09},
 \ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}):
 \begin{equation}
@@ -1151 +1172 @@
 (\ie fluxes plus content in mass exchanges).
 $\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter).
-Its default value is \np{rn\_atfp}~\forcode{= 10.e-3}.
+Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}.
 Note that the forcing correction term in the filter is not applied in linear free surface
-(\jp{lk\_vvl}~\forcode{= .false.}) (see \autoref{subsec:TRA_sbc}).
+(\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}).
 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$.
 
@@ -1166 +1187 @@
 % Equation of State (eosbn2) 
 % ================================================================
-\section{Equation of state (\protect\mdl{eosbn2}) }
+\section[Equation of state (\textit{eosbn2.F90})]
+{Equation of state (\protect\mdl{eosbn2})}
 \label{sec:TRA_eosbn2}
 %--------------------------------------------nameos-----------------------------------------------------
@@ -1176 +1198 @@
 %        Equation of State
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Equation of seawater (\protect\np{nn\_eos}~\forcode{= -1..1})}
+\subsection[Equation of seawater (\forcode{nn_eos = {-1,1}})]
+{Equation of seawater (\protect\np{nn\_eos}\forcode{ = {-1,1}})}
 \label{subsec:TRA_eos}
 
@@ -1186 +1209 @@
 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}.
+thus controlling rates of exchange between the atmosphere and the ocean interior \citep{roquet.madec.ea_JPO15}.
+Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{fofonoff.millard_bk83}) or
+TEOS-10 \citep{ioc.iapso_bk10} standards should be used anytime a simulation of the real ocean circulation is attempted
+\citep{roquet.madec.ea_JPO15}.
 The use of TEOS-10 is highly recommended because
 \textit{(i)}   it is the new official EOS,
@@ -1195 +1218 @@
 \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}.
+\citep{ioc.iapso_bk10, 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.
+To that purposed, a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is also available.
 
 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density.
@@ -1204 +1227 @@
 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 that value \citep{Gill1982}.
+density in the World Ocean varies by no more than 2$\%$ from that value \citep{gill_bk82}.
 
 Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which
@@ -1210 +1233 @@
 
 \begin{description}
-\item[\np{nn\_eos}~\forcode{= -1}]
-  the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used.
+\item[\np{nn\_eos}\forcode{ = -1}]
+  the polyTEOS10-bsq equation of seawater \citep{roquet.madec.ea_OM15} 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
@@ -1217 +1240 @@
   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}.
+  the TEOS-10 rational function approximation for hydrographic data analysis \citep{ioc.iapso_bk10}.
   A key point is that conservative state variables are used:
   Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{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}\,^{\circ}K^{-1}$, according to \citet{TEOS10}.
+  It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{ioc.iapso_bk10}.
   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
@@ -1229 +1252 @@
   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}~\forcode{= 0}]
+\item[\np{nn\_eos}\forcode{ = 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
@@ -1238 +1261 @@
   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}.
+  pressure \citep{fofonoff.millard_bk83}.
   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}~\forcode{= 1}]
-  a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen,
+\item[\np{nn\_eos}\forcode{ = 1}]
+  a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is chosen,
   the coefficients of which has been optimized to fit the behavior of TEOS10
-  (Roquet, personal comm.) (see also \citet{Roquet_JPO2015}).
+  (Roquet, personal comm.) (see also \citet{roquet.madec.ea_JPO15}).
   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}.
+  is enough for a proper treatment of the EOS in theoretical studies \citep{roquet.madec.ea_JPO15}.
   With such an equation of state there is no longer a distinction between
   \textit{conservative} and \textit{potential} temperature,
@@ -1303 +1326 @@
 %        Brunt-V\"{a}is\"{a}l\"{a} Frequency
 % -------------------------------------------------------------------------------------------------------------
-\subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}~\forcode{= 0..2})}
+\subsection[Brunt-V\"{a}is\"{a}l\"{a} frequency (\forcode{nn_eos = [0-2]})]
+{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}\forcode{ = [0-2]})}
 \label{subsec:TRA_bn2}
 
@@ -1329 +1353 @@
 \label{subsec:TRA_fzp}
 
-The freezing point of seawater is a function of salinity and pressure \citep{UNESCO1983}:
+The freezing point of seawater is a function of salinity and pressure \citep{fofonoff.millard_bk83}:
 \begin{equation}
   \label{eq:tra_eos_fzp}
@@ -1357 +1381 @@
 % Horizontal Derivative in zps-coordinate 
 % ================================================================
-\section{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})}
+\section[Horizontal derivative in \textit{zps}-coordinate (\textit{zpshde.F90})]
+{Horizontal derivative in \textit{zps}-coordinate (\protect\mdl{zpshde})}
 \label{sec:TRA_zpshde}
 
@@ -1363 +1388 @@
 I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln\_zps}~\forcode{= .true.}) at bottom and top (\np{ln\_isfcav}~\forcode{= .true.}),
+With partial cells (\np{ln\_zps}\forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{ = .true.}),
 in general, tracers in horizontally adjacent cells live at different depths.
 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and
 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module).
-The partial cell properties at the top (\np{ln\_isfcav}~\forcode{= .true.}) are computed in the same way as
+The partial cell properties at the top (\np{ln\_isfcav}\forcode{ = .true.}) are computed in the same way as
 for the bottom.
 So, only the bottom interpolation is explained below.
@@ -1379 +1404 @@
 \begin{figure}[!p]
   \begin{center}
-    \includegraphics[]{Fig_partial_step_scheme}
+    \includegraphics[width=\textwidth]{Fig_partial_step_scheme}
     \caption{
       \protect\label{fig:Partial_step_scheme}
       Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate
-      (\protect\np{ln\_zps}~\forcode{= .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$.
+      (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$.
       A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$,
       the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points.

@@ -1 @@
-*.aux
-*.bbl
-*.blg
-*.dvi
-*.fdb*
-*.fls
-*.idx
-*.ilg
-*.ind
-*.log
-*.maf
-*.mtc*
-*.out
-*.pdf
-*.toc
-_minted-*
+figures