Changeset 10502


Ignore:
Timestamp:
2019-01-10T18:45:21+01:00 (15 months ago)
Author:
nicolasmartin
Message:

Global work on math environnements for equations (partial commits)

Location:
NEMO/trunk/doc/latex/NEMO/subfiles
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

    r10442 r10502  
    33\begin{document} 
    44% ================================================================ 
    5 % Chapter 2 Space and Time Domain (DOM) 
     5% Chapter 2 ——— Space and Time Domain (DOM) 
    66% ================================================================ 
    77\chapter{Space Domain (DOM)} 
     
    4040\begin{figure}[!tb] 
    4141  \begin{center} 
    42     \includegraphics[width=0.90\textwidth]{Fig_cell} 
     42    \includegraphics[]{Fig_cell} 
    4343    \caption{ 
    4444      \protect\label{fig:cell} 
     
    4646      $t$ indicates scalar points where temperature, salinity, density, pressure and 
    4747      horizontal divergence are defined. 
    48       ($u$,$v$,$w$) indicates vector points, 
    49       and $f$ indicates vorticity points where both relative and planetary vorticities are defined 
     48      $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and 
     49      planetary vorticities are defined. 
    5050    } 
    5151  \end{center} 
     
    6464the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 
    6565 
    66 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 
    67 the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. 
     66The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 
     67gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    6868The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. 
    6969In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 
    7070the grid-point where the scale factors are defined. 
    7171Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 
    72 As a result, 
    73 the mesh on which partial derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, 
    74 and $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. 
    75 Discrete partial derivatives are formulated by the traditional, 
    76 centred second order finite difference approximation while 
    77 the scale factors are chosen equal to their local analytical value. 
     72As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
     73$\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. 
     74Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 
     75while the scale factors are chosen equal to their local analytical value. 
    7876An important point here is that the partial derivative of the scale factors must be evaluated by 
    7977centred finite difference approximation, not from their analytical expression. 
    80 This preserves the symmetry of the discrete set of equations and 
    81 therefore satisfies many of the continuous properties (see \autoref{apdx:C}). 
     78This preserves the symmetry of the discrete set of equations and therefore satisfies many of 
     79the continuous properties (see \autoref{apdx:C}). 
    8280A similar, related remark can be made about the domain size: 
    8381when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors 
    84 (see \autoref{eq:DOM_bar}) in the next section). 
     82(see \autoref{eq:DOM_bar} in the next section). 
    8583 
    8684%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    8987    \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 
    9088      \hline 
    91       T  &$i$     & $j$    & $k$     \\ \hline 
    92       u  & $i+1/2$   & $j$    & $k$    \\ \hline 
    93       v  & $i$    & $j+1/2$   & $k$    \\ \hline 
    94       w  & $i$    & $j$    & $k+1/2$   \\ \hline 
    95       f  & $i+1/2$   & $j+1/2$   & $k$    \\ \hline 
    96       uw & $i+1/2$   & $j$    & $k+1/2$   \\ \hline 
    97       vw & $i$    & $j+1/2$   & $k+1/2$   \\ \hline 
    98       fw & $i+1/2$   & $j+1/2$   & $k+1/2$   \\ \hline 
     89      T  & $i      $ & $j      $ & $k      $ \\ 
     90      \hline 
     91      u  & $i + 1/2$ & $j      $ & $k      $ \\ 
     92      \hline 
     93      v  & $i      $ & $j + 1/2$ & $k      $ \\ 
     94      \hline 
     95      w  & $i      $ & $j      $ & $k + 1/2$ \\ 
     96      \hline 
     97      f  & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
     98      \hline 
     99      uw & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
     100      \hline 
     101      vw & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
     102      \hline 
     103      fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
     104      \hline 
    99105    \end{tabular} 
    100106    \caption{ 
    101107      \protect\label{tab:cell} 
    102108      Location of grid-points as a function of integer or integer and a half value of the column, line or level. 
    103       This indexing is only used for the writing of the semi-discrete equation. 
     109      This indexing is only used for the writing of the semi -discrete equation. 
    104110      In the code, the indexing uses integer values only and has a reverse direction in the vertical 
    105111      (see \autoref{subsec:DOM_Num_Index}) 
     
    115121\label{subsec:DOM_operators} 
    116122 
    117 Given the values of a variable $q$ at adjacent points, 
    118 the differencing and averaging operators at the midpoint between them are: 
    119 \[ 
     123Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at 
     124the midpoint between them are: 
     125\begin{alignat*}{2} 
    120126  % \label{eq:di_mi} 
    121   \begin{split} 
    122     \delta_i [q]       &=  \  \    q(i+1/2)  - q(i-1/2)     \\ 
    123     \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2 
    124   \end{split} 
    125 \] 
    126  
    127 Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$. 
     127  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\ 
     128  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 
     129\end{alignat*} 
     130 
     131Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 
    128132Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 
    129133a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 
    130 its Laplacien is defined at $t$-point. 
    131 These operators have the following discrete forms in the curvilinear $s$-coordinate system: 
     134its Laplacian is defined at $t$-point. 
     135These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
    132136\[ 
    133137  % \label{eq:DOM_grad} 
    134   \nabla q\equiv  \frac{1}{e_{1u} } \delta_{i+1/2 } [q] \;\,\mathbf{i} 
    135   +   \frac{1}{e_{2v} } \delta_{j+1/2 } [q] \;\,\mathbf{j} 
    136   +   \frac{1}{e_{3w}} \delta_{k+1/2} [q] \;\,\mathbf{k} 
     138  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i 
     139                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j 
     140                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 
    137141\] 
    138142\begin{multline*} 
    139143  % \label{eq:DOM_lap} 
    140   \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    141   \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] 
    142     +                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)     \\ 
    143   +\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right] 
     144  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
     145                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) 
     146                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ 
     147                  + \frac{1}{e_{3t}} 
     148                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 
    144149\end{multline*} 
    145150 
    146 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ 
    147 defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, 
    148 and its divergence defined at $t$-points: 
    149 \begin{align*} 
    150   % \label{eq:DOM_curl} 
    151   \nabla \times {\rm{\bf A}}\equiv & 
    152                                      \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \mathbf{i} \\ 
    153   +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \mathbf{j} \\ 
    154   +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \mathbf{k} 
    155 \end{align*} 
    156 \begin{align*} 
    157   % \label{eq:DOM_div} 
    158   \nabla \cdot \rm{\bf A} \equiv 
    159   \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 
    160   +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 
    161 \end{align*} 
     151Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
     152vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 
     153its divergence defined at $t$-points: 
     154\begin{multline} 
     155% \label{eq:DOM_curl} 
     156  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} 
     157                                 \Big[   \delta_{j + 1/2} (e_{3w} \, a_3) 
     158                                       - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\ 
     159                               + \frac{1}{e_{2u} \, e_{3uw}} 
     160                                 \Big[   \delta_{k + 1/2} (e_{1u} \, a_1) 
     161                                       - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\ 
     162                               + \frac{1}{e_{1f} \, e_{2f}} 
     163                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) 
     164                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 
     165\end{multline} 
     166\begin{equation} 
     167% \label{eq:DOM_div} 
     168  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
     169                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 
     170                              + \frac{1}{e_{3t}} \delta_k (a_3) 
     171\end{equation} 
    162172 
    163173The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which 
    164 is a masked field (\ie equal to zero inside solid area): 
     174is a masked field (i.e. equal to zero inside solid area): 
    165175\begin{equation} 
    166176  \label{eq:DOM_bar} 
    167   \bar q    =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk} 
    168   \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } 
     177  \bar q = \frac{1}{H} \int_{k^b}^{k^o} q \; e_{3q} \, dk \equiv \frac{1}{H_q} \sum \limits_k q \; e_{3q} 
    169178\end{equation} 
    170179where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 
    171 $k^b$ and $k^o$ are the bottom and surface $k$-indices, 
    172 and the symbol $k^o$ refers to a summation over all grid points of the same type in the direction indicated by 
    173 the subscript (here $k$).  
     180$k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over 
     181all grid points of the same type in the direction indicated by the subscript (here $k$). 
    174182 
    175183In continuous form, the following properties are satisfied: 
    176 \begin{equation} 
     184\begin{gather} 
    177185  \label{eq:DOM_curl_grad} 
    178   \nabla \times \nabla q ={\rm {\bf {0}}} 
    179 \end{equation} 
    180 \begin{equation} 
     186  \nabla \times \nabla q = \vect 0 \\ 
    181187  \label{eq:DOM_div_curl} 
    182   \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0 
    183 \end{equation} 
     188  \nabla \cdot (\nabla \times \vect A) = 0 
     189\end{gather} 
    184190 
    185191It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as 
    186 the scalar $q$ is taken at $t$-points and 
    187 the vector \textbf{A} has its components defined at vector points $(u,v,w)$. 
     192the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 
     193vector points $(u,v,w)$. 
    188194 
    189195Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. 
    190 Using integration by parts it can be shown that 
    191 the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, 
    192 and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, $\overline{\,\cdot\,}^{\,k}$ and 
    193 $\overline{\,\cdot\,}^{\,k}$) are symmetric linear operators, \ie 
    194 \begin{align} 
     196Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
     197are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 
     198$\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
     199\begin{alignat}{4} 
    195200  \label{eq:DOM_di_adj} 
    196   \sum\limits_i { a_i \;\delta_i \left[ b \right]} 
    197   &\equiv -\sum\limits_i {\delta_{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\ 
     201  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\ 
    198202  \label{eq:DOM_mi_adj} 
    199   \sum\limits_i { a_i \;\overline b^{\,i}} 
    200   & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 
    201 \end{align} 
    202  
    203 In other words, the adjoint of the differencing and averaging operators are $\delta_i^*=\delta_{i+1/2}$ and  
    204 ${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively.  
     203  &\sum \limits_i a_i \; \overline b^{\, i} &\equiv &  &&\sum \limits_i \overline a ^{\, i + 1/2}     &b_{i + 1/2} 
     204\end{alignat} 
     205 
     206In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and  
     207$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 
    205208These two properties will be used extensively in the \autoref{apdx:C} to 
    206209demonstrate integral conservative properties of the discrete formulation chosen. 
     
    215218\begin{figure}[!tb] 
    216219  \begin{center} 
    217     \includegraphics[width=0.90\textwidth]{Fig_index_hor} 
     220    \includegraphics[]{Fig_index_hor} 
    218221    \caption{ 
    219222      \protect\label{fig:index_hor} 
     
    230233Therefore a specific integer indexing must be defined for points other than $t$-points 
    231234(\ie velocity and vorticity grid-points). 
    232 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k=1$. 
     235Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. 
    233236 
    234237% ----------------------------------- 
     
    250253\label{subsec:DOM_Num_Index_vertical} 
    251254 
    252 In the vertical, the chosen indexing requires special attention since 
    253 the $k$-axis is re-orientated downward in the \fortran code compared to 
    254 the indexing used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 
    255 The sea surface corresponds to the $w$-level $k=1$ which is the same index as $t$-level just below 
     255In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in 
     256the \fortran code compared to the indexing used in the semi -discrete equations and 
     257given in \autoref{subsec:DOM_cell}. 
     258The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below 
    256259(\autoref{fig:index_vert}). 
    257 The last $w$-level ($k=jpk$) either corresponds to the ocean floor or is inside the bathymetry while 
     260The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 
    258261the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 
    259262Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, 
     
    262265have the same $i$ or $j$ index 
    263266(compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 
    264 Since the scale factors are chosen to be strictly positive, a \emph{minus sign} appears in the \fortran  
    265 code \emph{before all the vertical derivatives} of the discrete equations given in this documentation. 
     267Since the scale factors are chosen to be strictly positive, 
     268a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of 
     269the discrete equations given in this documentation. 
    266270 
    267271%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    268272\begin{figure}[!pt] 
    269273  \begin{center} 
    270     \includegraphics[width=.90\textwidth]{Fig_index_vert} 
     274    \includegraphics[]{Fig_index_vert} 
    271275    \caption{ 
    272276      \protect\label{fig:index_vert} 
     
    287291The total size of the computational domain is set by the parameters \np{jpiglo}, 
    288292\np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 
    289 %%% 
    290 %%% 
    291 %%% 
    292293Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 
    293294the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 
     
    299300\section{Needed fields} 
    300301\label{sec:DOM_fields} 
    301 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
     302The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 
     303gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    302304The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 
    303305The associated scale factors are defined using the analytical first derivative of the transformation 
    304306\autoref{eq:scale_factors}. 
    305 Necessary fields for configuration definition are: \\ 
    306 Geographic position : 
    307  
    308 longitude: glamt, glamu, glamv and glamf (at T, U, V and F point) 
    309  
    310 latitude: gphit, gphiu, gphiv and gphif (at T, U, V and F point)\\ 
    311 Coriolis parameter (if domain not on the sphere):  
    312  
    313  ff\_f  and  ff\_t (at T and F point)\\ 
    314 Scale factors :  
     307Necessary fields for configuration definition are: 
     308 
     309\begin{itemize} 
     310\item 
     311  Geographic position: 
     312  longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 
     313  latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 
     314  (all respectively at T, U, V and F point) 
     315\item 
     316  Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 
     317  (at T and F point) 
     318\item 
     319  Scale factors: 
     320  \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 
     321  \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 
     322  \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 
     323  \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits),  
     324  \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 
     325\end{itemize} 
    315326  
    316  e1t, e1u, e1v and e1f (on i direction), 
    317  
    318  e2t, e2u, e2v and e2f (on j direction) and 
    319  
    320  ie1e2u\_v, e1e2u , e1e2v    
    321   
    322 e1e2u , e1e2v are u and v surfaces (if gridsize reduction in some straits)\\ 
    323 ie1e2u\_v is a flag to flag set u and  v surfaces are neither read nor computed.\\ 
    324   
    325 These fields can be read in an domain input file which name is setted in 
    326 \np{cn\_domcfg} parameter specified in \ngn{namcfg}. 
     327These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 
     328\ngn{namcfg}. 
    327329 
    328330\nlst{namcfg} 
    329 or they can be defined in an analytical way in MY\_SRC directory of the configuration. 
     331 
     332Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 
    330333For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 
    331 For analytical definition of input fields two routines are supplied: \mdl{userdef\_hgr} and \mdl{userdef\_zgr}. 
    332 They are an example of GYRE configuration parameters, and they are available in NEMO/OPA\_SRC/USR directory, 
    333 they provide the horizontal and vertical mesh.  
     334For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
     335They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 
     336they provide the horizontal and vertical mesh. 
    334337% ------------------------------------------------------------------------------------------------------------- 
    335338%        Needed fields  
     
    366369($i$ and $j$, respectively) (geographical configuration of the mesh), 
    367370the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 
    368 and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the \mdl{domhgr} module. 
     371and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 
    369372The model computes the grid-point positions and scale factors in the horizontal plane as follows: 
    370 \begin{flalign*} 
    371   \lambda_t &\equiv \text{glamt}= \lambda(i)   & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ 
    372   \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ 
    373   \lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ 
    374   \lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) 
    375 \end{flalign*} 
    376 \begin{flalign*} 
    377   e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)     \; \cos\varphi(j)  |& 
    378   e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ 
    379   e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2) \; \cos\varphi(j)  |& 
    380   e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ 
    381   e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i)     \; \cos\varphi(j+1/2)  |& 
    382   e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\ 
    383   e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2)  |& 
    384   e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)| 
    385 \end{flalign*} 
     373\begin{align*} 
     374   \lambda_t &\equiv \text{glamt} =      \lambda (i      ) 
     375  &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\ 
     376   \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2) 
     377  &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\ 
     378   \lambda_v &\equiv \text{glamv} =      \lambda (i      ) 
     379  &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\ 
     380   \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2) 
     381  &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\ 
     382   e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) | 
     383  &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
     384   e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) | 
     385  &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
     386   e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) | 
     387  &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\ 
     388   e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 
     389  &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | 
     390\end{align*} 
    386391where the last letter of each computational name indicates the grid point considered and 
    387392$r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 
    388393Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 
    389 thus no specific arrays are defined at $w$-points.  
     394thus no specific arrays are defined at $w$-points. 
    390395 
    391396Note that the definition of the scale factors 
     
    405410\begin{figure}[!t] 
    406411  \begin{center} 
    407     \includegraphics[width=0.90\textwidth]{Fig_zgr_e3} 
     412    \includegraphics[]{Fig_zgr_e3} 
    408413    \caption{ 
    409414      \protect\label{fig:zgr_e3} 
    410415      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
    411416      and (b) analytically derived grid-point position and scale factors. 
    412       For both grids here, 
    413       the same $w$-point depth has been chosen but in (a) the $t$-points are set half way between $w$-points while 
    414       in (b) they are defined from an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$. 
     417      For both grids here, the same $w$-point depth has been chosen but 
     418      in (a) the $t$-points are set half way between $w$-points while 
     419      in (b) they are defined from an analytical function: 
     420      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 
    415421      Note the resulting difference between the value of the grid-size $\Delta_k$ and 
    416422      those of the scale factor $e_k$. 
     
    426432\label{subsec:DOM_hgr_msh_choice} 
    427433 
    428  
    429434% ------------------------------------------------------------------------------------------------------------- 
    430435%        Grid files 
     
    434439 
    435440All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 
    436 can be saved in files if \np{nn\_msh} $\not= 0$ (namelist variable in \ngn{namdom}). 
     441can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 
    437442This can be particularly useful for plots and off-line diagnostics. 
    438443In some cases, the user may choose to make a local modification of a scale factor in the code. 
     
    441446An example is Gibraltar Strait in the ORCA2 configuration. 
    442447When such modifications are done, 
    443 the output grid written when \np{nn\_msh} $\not= 0$ is no more equal to the input grid. 
     448the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 
    444449 
    445450% ================================================================ 
     
    466471\begin{figure}[!tb] 
    467472  \begin{center} 
    468     \includegraphics[width=1.0\textwidth]{Fig_z_zps_s_sps} 
     473    \includegraphics[]{Fig_z_zps_s_sps} 
    469474    \caption{ 
    470475      \protect\label{fig:z_zps_s_sps} 
     
    475480      (d) hybrid $s-z$ coordinate, 
    476481      (e) hybrid $s-z$ coordinate with partial step, and 
    477       (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}). 
     482      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). 
    478483      Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 
    479484    } 
     
    485490must be done once of all at the beginning of an experiment. 
    486491It is not intended as an option which can be enabled or disabled in the middle of an experiment. 
    487 Three main choices are offered (\autoref{fig:z_zps_s_sps}a to c): 
    488 $z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}), 
    489 $z$-coordinate with partial step bathymetry (\np{ln\_zps}\forcode{ = .true.}), 
    490 or generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}). 
     492Three main choices are offered (\autoref{fig:z_zps_s_sps}): 
     493$z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 
     494$z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 
     495or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 
    491496Hybridation of the three main coordinates are available: 
    492 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}e). 
     497$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 
    493498By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 
    494 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}f). 
    495 When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}), 
    496 the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface 
    497 (in other words, the top of the ocean in not a rigid-lid).  
     499the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 
     500When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 
     501the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 
     502(in other words, the top of the ocean in not a rigid-lid). 
    498503The last choice in terms of vertical coordinate concerns the presence (or not) in 
    499504the model domain of ocean cavities beneath ice shelves. 
     
    502507and partial step are also applied at the ocean/ice shelf interface. 
    503508 
    504 Contrary to the horizontal grid, the vertical grid is computed in the code and 
    505 no provision is made for reading it from a file. 
     509Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 
     510reading it from a file. 
    506511The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 
    507512\footnote{ 
    508513  N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 
    509   so that the computation of the number of wet ocean point in each water column is by-passed 
    510 }.  
    511 If \np{ln\_isfcav}\forcode{ = .true.}, 
    512 an extra file input file describing the ice shelf draft (in meters) (\ifile{isf\_draft\_meter}) is needed. 
     514  so that the computation of the number of wet ocean point in each water column is by-passed}. 
     515If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 
     516the ice shelf draft (in meters) is needed. 
    513517 
    514518After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 
    515519\begin{description} 
    516520\item[\textit{zco}] 
    517   set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$. 
     521  set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 
    518522\item[\textit{zps}] 
    519   set a reference coordinate transformation $z_0 (k)$, 
    520   and calculate the thickness of the deepest level at each $(i,j)$ point using the bathymetry, 
    521   to obtain the final three-dimensional depth and scale factor arrays. 
     523  set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 
     524  each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 
    522525\item[\textit{sco}] 
    523   smooth the bathymetry to fulfil the hydrostatic consistency criteria and 
     526  smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    524527  set the three-dimensional transformation. 
    525528\item[\textit{s-z} and \textit{s-zps}] 
    526   smooth the bathymetry to fulfil the hydrostatic consistency criteria and 
     529  smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    527530  set the three-dimensional transformation $z(i,j,k)$, 
    528531  and possibly introduce masking of extra land points to better fit the original bathymetry file. 
     
    532535%%% 
    533536 
    534 Unless a linear free surface is used (\np{ln\_linssh}\forcode{ = .false.}), 
     537Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 
    535538the arrays describing the grid point depths and vertical scale factors are three set of 
    536539three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 
    537 The time at which they are defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively. 
     540The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
    538541They are updated at each model time step using a fixed reference coordinate system which 
    539542computer names have a $\_0$ suffix. 
    540 When the linear free surface option is used (\np{ln\_linssh}\forcode{ = .true.}), 
    541 \textit{before}, \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart.  
    542  
     543When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 
     544\textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 
    543545 
    544546% ------------------------------------------------------------------------------------------------------------- 
     
    551553(found in \ngn{namdom} namelist):  
    552554\begin{description} 
    553 \item[\np{nn\_bathy}\forcode{ = 0}]: 
     555\item[\np{nn\_bathy}~\forcode{= 0}]: 
    554556  a flat-bottom domain is defined. 
    555557  The total depth $z_w (jpk)$ is given by the coordinate transformation. 
    556   The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}.  
    557 \item[\np{nn\_bathy}\forcode{ = -1}]: 
     558  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 
     559\item[\np{nn\_bathy}~\forcode{= -1}]: 
    558560  a domain with a bump of topography one third of the domain width at the central latitude. 
    559   This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount.  
    560 \item[\np{nn\_bathy}\forcode{ = 1}]: 
     561  This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 
     562\item[\np{nn\_bathy}~\forcode{= 1}]: 
    561563  read a bathymetry and ice shelf draft (if needed). 
    562564  The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at 
     
    569571  The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 
    570572  each grid point of the model grid. 
    571   This file is only needed if \np{ln\_isfcav}\forcode{ = .true.}. 
     573  This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 
    572574  Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 
    573575\end{description} 
    574576 
    575577When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 
    576 (e.g, great lakes, Caspian sea...) 
    577 even if the model resolution does not allow their communication with the rest of the ocean. 
     578(\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 
     579the rest of the ocean. 
    578580This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
    579581so these seas can be removed from the ocean domain. 
    580582The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 
    581 but the code has to be adapted to the user's configuration.  
     583but the code has to be adapted to the user's configuration. 
    582584 
    583585% ------------------------------------------------------------------------------------------------------------- 
    584586%        z-coordinate  and reference coordinate transformation 
    585587% ------------------------------------------------------------------------------------------------------------- 
    586 \subsection[$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and ref. coordinate] 
    587             {$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and reference coordinate} 
     588\subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 
     589            {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 
    588590\label{subsec:DOM_zco} 
    589591 
     
    591593\begin{figure}[!tb] 
    592594  \begin{center} 
    593     \includegraphics[width=0.90\textwidth]{Fig_zgr} 
     595    \includegraphics[]{Fig_zgr} 
    594596    \caption{ 
    595597      \protect\label{fig:zgr} 
     
    602604%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    603605 
    604 The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ and $gdepw_0$ for 
    605 $t$- and $w$-points, respectively. 
    606 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the ocean surface. 
     606The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 
     607respectively. 
     608As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 
     609$gdepw_0(1)$ is the ocean surface. 
    607610There are at most \jp{jpk}-1 $t$-points inside the ocean, 
    608 the additional $t$-point at $jk=jpk$ is below the sea floor and is not used. 
     611the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 
    609612The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 
    610613analytical derivative with respect to $k$ provides the vertical scale factors. 
     
    613616using parameters provided in the \ngn{namcfg} namelist. 
    614617 
    615 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr=0}). 
    616 In that case, 
    617 the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} (total ocean depth in meters) fully define the grid.  
     618It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 
     619In that case, the parameters \jp{jpk} (number of $w$-levels) and 
     620\np{pphmax} (total ocean depth in meters) fully define the grid. 
    618621 
    619622For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 
    620623The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps):  
    621 \begin{equation} 
     624\begin{gather} 
    622625  \label{eq:DOM_zgr_ana_1} 
    623   \begin{split} 
    624     z_0 (k)    &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 
    625     e_3^0 (k)  &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right| 
    626   \end{split} 
    627 \end{equation} 
    628 where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. 
     626    z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 
     627    e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt| 
     628\end{gather} 
     629where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 
    629630Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 
    630631a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 
    631632 
    632 If the ice shelf cavities are opened (\np{ln\_isfcav}\forcode{ = .true.}), the definition of $z_0$ is the same. 
     633If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 
    633634However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 
    634635\begin{equation} 
    635636  \label{eq:DOM_zgr_ana_2} 
    636637  \begin{split} 
    637     e_3^T(k) &= z_W (k+1) - z_W (k)  \\ 
    638     e_3^W(k) &= z_T (k)   - z_T (k-1) \\ 
     638    e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\ 
     639    e_3^W(k) &= z_T (k    ) - z_T (k - 1) 
    639640  \end{split} 
    640641\end{equation} 
    641642This formulation decrease the self-generated circulation into the ice shelf cavity  
    642643(which can, in extreme case, leads to blow up).\\ 
    643  
    644644  
    645 The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the surface (bottom) layers and 
     645The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 
    646646a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 
    647647This leads to the following conditions: 
    648648\begin{equation} 
    649649  \label{eq:DOM_zgr_coef} 
    650   \begin{split} 
    651     e_3 (1+1/2)      &=10. \\ 
    652     e_3 (jpk-1/2) &=500. \\ 
    653     z(1)       &=0. \\ 
    654     z(jpk)        &=-5000. \\ 
    655   \end{split} 
     650  \begin{array}{ll} 
     651    e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\ 
     652    e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 
     653  \end{array} 
    656654\end{equation} 
    657655 
    658 With the choice of the stretching $h_{cr} =3$ and the number of levels \jp{jpk}=$31$, 
    659 the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in 
     656With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 
     657the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 
    660658\autoref{eq:DOM_zgr_ana_2} have been determined such that 
    661659\autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 
    662660For the first standard ORCA2 vertical grid this led to the following values: 
    663 $h_{sur} =4762.96$, $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. 
     661$h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 
    664662The resulting depths and scale factors as a function of the model levels are shown in 
    665663\autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 
    666 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.  
    667  
    668 Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is possible to recalculate them. 
    669 In that case the user sets \np{ppsur}\forcode{ = }\np{ppa0}\forcode{ = }\np{ppa1}\forcode{ = 999999}., 
     664Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 
     665 
     666Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 
     667In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 
    670668in \ngn{namcfg} namelist, and specifies instead the four following parameters: 
    671669\begin{itemize} 
    672670\item 
    673   \np{ppacr}=$h_{cr} $: stretching factor (nondimensional). 
     671  \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 
    674672  The larger \np{ppacr}, the smaller the stretching. 
    675673  Values from $3$ to $10$ are usual. 
    676674\item 
    677   \np{ppkth}=$h_{th} $: is approximately the model level at which maximum stretching occurs 
     675  \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 
    678676  (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 
    679677\item 
     
    683681\end{itemize} 
    684682As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 
    685 \jp{jpk}\forcode{ = 46}, \np{ppacr}\forcode{ = 9}, \np{ppkth}\forcode{ = 23.563}, 
    686 \np{ppdzmin}\forcode{ = 6}m, \np{pphmax}\forcode{ = 5750}m. 
     683\jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 
    687684 
    688685%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    691688    \begin{tabular}{c||r|r|r|r} 
    692689      \hline 
    693       \textbf{LEVEL}& \textbf{gdept\_1d}& \textbf{gdepw\_1d}& \textbf{e3t\_1d }& \textbf{e3w\_1d  } \\ \hline 
    694       1  &  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline 
    695       2  &  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline 
    696       3  &  \textbf{25.00} &    20.00 &   \textbf{ 10.00} &     10.00 \\   \hline 
    697       4  &  \textbf{35.01} &    30.00 &   \textbf{ 10.01} &     10.00 \\   \hline 
    698       5  &  \textbf{45.01} &    40.01 &   \textbf{ 10.01} &  10.01 \\   \hline 
    699       6  &  \textbf{55.03} &    50.02 &   \textbf{ 10.02} &     10.02 \\   \hline 
    700       7  &  \textbf{65.06} &    60.04 &   \textbf{ 10.04} &  10.03 \\   \hline 
    701       8  &  \textbf{75.13} &    70.09 &   \textbf{ 10.09} &  10.06 \\   \hline 
    702       9  &  \textbf{85.25} &    80.18 &   \textbf{ 10.17} &  10.12 \\   \hline 
    703       10 &  \textbf{95.49} &    90.35 &   \textbf{ 10.33} &  10.24 \\   \hline 
    704       11 &  \textbf{105.97}   &   100.69 &   \textbf{ 10.65} &  10.47 \\   \hline 
    705       12 &  \textbf{116.90}   &   111.36 &   \textbf{ 11.27} &  10.91 \\   \hline 
    706       13 &  \textbf{128.70}   &   122.65 &   \textbf{ 12.47} &  11.77 \\   \hline 
    707       14 &  \textbf{142.20}   &   135.16 &   \textbf{ 14.78} &  13.43 \\   \hline 
    708       15 &  \textbf{158.96}   &   150.03 &   \textbf{ 19.23} &  16.65 \\   \hline 
    709       16 &  \textbf{181.96}   &   169.42 &   \textbf{ 27.66} &  22.78 \\   \hline 
    710       17 &  \textbf{216.65}   &   197.37 &   \textbf{ 43.26} &  34.30 \\ \hline 
    711       18 &  \textbf{272.48}   &   241.13 &   \textbf{ 70.88} &  55.21 \\ \hline 
    712       19 &  \textbf{364.30}   &   312.74 &   \textbf{116.11} &  90.99 \\ \hline 
    713       20 &  \textbf{511.53}   &   429.72 &   \textbf{181.55} &    146.43 \\ \hline 
    714       21 &  \textbf{732.20}   &   611.89 &   \textbf{261.03} &    220.35 \\ \hline 
    715       22 &  \textbf{1033.22}&  872.87 &   \textbf{339.39} &    301.42 \\ \hline 
    716       23 &  \textbf{1405.70}& 1211.59 & \textbf{402.26} &   373.31 \\ \hline 
    717       24 &  \textbf{1830.89}& 1612.98 & \textbf{444.87} &   426.00 \\ \hline 
    718       25 &  \textbf{2289.77}& 2057.13 & \textbf{470.55} &   459.47 \\ \hline 
    719       26 &  \textbf{2768.24}& 2527.22 & \textbf{484.95} &   478.83 \\ \hline 
    720       27 &  \textbf{3257.48}& 3011.90 & \textbf{492.70} &   489.44 \\ \hline 
    721       28 &  \textbf{3752.44}& 3504.46 & \textbf{496.78} &   495.07 \\ \hline 
    722       29 &  \textbf{4250.40}& 4001.16 & \textbf{498.90} &   498.02 \\ \hline 
    723       30 &  \textbf{4749.91}& 4500.02 & \textbf{500.00} &   499.54 \\ \hline 
    724       31 &  \textbf{5250.23}& 5000.00 &   \textbf{500.56} & 500.33 \\ \hline 
     690      \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 
     691      \hline 
     692      1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\ 
     693      \hline 
     694      2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\ 
     695      \hline 
     696      3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\ 
     697      \hline 
     698      4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\ 
     699      \hline 
     700      5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\ 
     701      \hline 
     702      6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\ 
     703      \hline 
     704      7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\ 
     705      \hline 
     706      8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\ 
     707      \hline 
     708      9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\ 
     709      \hline 
     710      10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\ 
     711      \hline 
     712      11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\ 
     713      \hline 
     714      12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\ 
     715      \hline 
     716      13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\ 
     717      \hline 
     718      14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\ 
     719      \hline 
     720      15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\ 
     721      \hline 
     722      16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\ 
     723      \hline 
     724      17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\ 
     725      \hline 
     726      18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\ 
     727      \hline 
     728      19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\ 
     729      \hline 
     730      20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\ 
     731      \hline 
     732      21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\ 
     733      \hline 
     734      22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\ 
     735      \hline 
     736      23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\ 
     737      \hline 
     738      24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\ 
     739      \hline 
     740      25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\ 
     741      \hline 
     742      26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\ 
     743      \hline 
     744      27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\ 
     745      \hline 
     746      28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\ 
     747      \hline 
     748      29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\ 
     749      \hline 
     750      30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\ 
     751      \hline 
     752      31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\ 
     753      \hline 
    725754    \end{tabular} 
    726755  \end{center} 
     
    736765%        z-coordinate with partial step 
    737766% ------------------------------------------------------------------------------------------------------------- 
    738 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}\forcode{ = .true.})} 
     767\subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 
    739768\label{subsec:DOM_zps} 
    740769%--------------------------------------------namdom------------------------------------------------------- 
     
    744773 
    745774In $z$-coordinate partial step, 
    746 the depths of the model levels are defined by the reference analytical function $z_0 (k)$ as described in 
    747 the previous section, \emph{except} in the bottom layer. 
     775the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 
     776the previous section, \textit{except} in the bottom layer. 
    748777The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 
    749778allow a better representation of the bathymetry, especially in the case of small slopes 
     
    752781With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 
    753782The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 
    754 the maximum thickness allowed is $2*e_{3t}(jpk-1)$. 
     783the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 
    755784This has to be kept in mind when specifying values in \ngn{namdom} namelist, 
    756785as the maximum depth \np{pphmax} in partial steps: 
    757 for example, with \np{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, 
    758 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). 
     786for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 
     787the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 
    759788Two variables in the namdom namelist are used to define the partial step vertical grid. 
    760789The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is 
     
    767796%        s-coordinate 
    768797% ------------------------------------------------------------------------------------------------------------- 
    769 \subsection{$S$-coordinate (\protect\np{ln\_sco}\forcode{ = .true.})} 
     798\subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 
    770799\label{subsec:DOM_sco} 
    771800%------------------------------------------nam_zgr_sco--------------------------------------------------- 
     
    774803%-------------------------------------------------------------------------------------------------------------- 
    775804Options are defined in \ngn{namzgr\_sco}. 
    776 In $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}), the depth and thickness of the model levels are defined from 
     805In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 
    777806the product of a depth field and either a stretching function or its derivative, respectively: 
    778807 
    779 \[ 
     808\begin{align*} 
    780809  % \label{eq:DOM_sco_ana} 
    781   \begin{split} 
    782     z(k)       &= h(i,j) \; z_0(k)  \\ 
    783     e_3(k)  &= h(i,j) \; z_0'(k) 
    784   \end{split} 
    785 \] 
     810  z(k)   &= h(i,j) \; z_0 (k) \\ 
     811  e_3(k) &= h(i,j) \; z_0'(k) 
     812\end{align*} 
    786813 
    787814where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and 
     
    789816The depth field $h$ is not necessary the ocean depth, 
    790817since a mixed step-like and bottom-following representation of the topography can be used 
    791 (\autoref{fig:z_zps_s_sps}d-e) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}f). 
     818(\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 
    792819The namelist parameter \np{rn\_rmax} determines the slope at which 
    793 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.  
     820the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 
    794821The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 
    795822the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 
     
    799826 
    800827The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 
    801 (\np{ln\_s\_SH94}\forcode{ = .false.} and \np{ln\_s\_SF12}\forcode{ = .false.}).  
     828(\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 
    802829This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 
    803830 
    804831\[ 
    805   z = s_{min}+C\left(s\right)\left(H-s_{min}\right) 
     832  z = s_{min} + C (s) (H - s_{min}) 
    806833  % \label{eq:SH94_1} 
    807834\] 
     
    810837allows a $z$-coordinate to placed on top of the stretched coordinate, 
    811838and $z$ is the depth (negative down from the asea surface). 
     839\begin{gather*} 
     840  s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 
     841  % \label{eq:DOM_s} 
     842 \\ 
     843  % \label{eq:DOM_sco_function} 
     844  C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 
     845\end{gather*} 
     846 
     847A stretching function, 
     848modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 
     849is also available and is more commonly used for shelf seas modelling: 
    812850 
    813851\[ 
    814   s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1 
    815   % \label{eq:DOM_s} 
    816 \] 
    817  
    818 \[ 
    819   % \label{eq:DOM_sco_function} 
    820   \begin{split} 
    821     C(s) &=  \frac{ \left[   \tanh{ \left( \theta \, (s+b) \right)} 
    822         - \tanh{ \left(  \theta \, b      \right)}  \right]} 
    823     {2\;\sinh \left( \theta \right)} 
    824   \end{split} 
    825 \] 
    826  
    827 A stretching function, 
    828 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}\forcode{ = .true.}), 
    829 is also available and is more commonly used for shelf seas modelling: 
    830  
    831 \[ 
    832   C\left(s\right) =   \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} +      \\ 
    833   b\frac{ \tanh \left[ \theta \left(s + \frac{1}{2} \right)\right] - \tanh\left(\frac{\theta}{2}\right)}{ 2\tanh\left (\frac{\theta}{2}\right)} 
     852  C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 
     853         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)} 
     854                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)} 
    834855  % \label{eq:SH94_2} 
    835856\] 
     
    838859\begin{figure}[!ht] 
    839860  \begin{center} 
    840     \includegraphics[width=1.0\textwidth]{Fig_sco_function} 
     861    \includegraphics[]{Fig_sco_function} 
    841862    \caption{ 
    842863      \protect\label{fig:sco_function} 
     
    848869%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    849870 
    850 where $H_c$ is the critical depth (\np{rn\_hc}) at which 
    851 the coordinate transitions from pure $\sigma$ to the stretched coordinate, 
    852 and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and bottom control parameters such that 
    853 $0\leqslant \theta \leqslant 20$, and $0\leqslant b\leqslant 1$. 
     871where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 
     872the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 
     873bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 
    854874$b$ has been designed to allow surface and/or bottom increase of the vertical resolution 
    855875(\autoref{fig:sco_function}). 
     
    859879In this case the a stretching function $\gamma$ is defined such that: 
    860880 
    861 \[ 
    862   z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1 
     881\begin{equation} 
     882  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 
    863883  % \label{eq:z} 
    864 \] 
     884\end{equation} 
    865885 
    866886The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 
    867887 
    868 \[ 
     888\begin{gather*} 
    869889  % \label{eq:DOM_gamma_deriv} 
    870   \gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right) 
    871 \] 
    872  
    873 Where: 
    874 \[ 
     890  \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt) 
     891           + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\ 
     892  \intertext{Where:} 
    875893  % \label{eq:DOM_gamma} 
    876   f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1} 
    877 \] 
     894  f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 
     895  \quad \text{and} \quad \sigma = \frac{k}{n - 1} 
     896\end{gather*} 
    878897 
    879898This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 
     
    892911%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    893912\begin{figure}[!ht] 
    894    \includegraphics[width=1.0\textwidth]{Fig_DOM_compare_coordinates_surface} 
    895    \caption{ 
    896      A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 
    897      a 50 level $Z$-coordinate (contoured surfaces) and 
    898      the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in 
    899      the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. 
    900      For clarity every third coordinate surface is shown. 
    901    } 
    902    \label{fig:fig_compare_coordinates_surface} 
     913  \includegraphics[]{Fig_DOM_compare_coordinates_surface} 
     914  \caption{ 
     915    A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 
     916    a 50 level $Z$-coordinate (contoured surfaces) and 
     917    the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for 
     918    a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 
     919    For clarity every third coordinate surface is shown. 
     920  } 
     921  \label{fig:fig_compare_coordinates_surface} 
    903922\end{figure} 
    904 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     923 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    905924 
    906925This gives a smooth analytical stretching in computational space that is constrained to 
     
    925944 
    926945% ------------------------------------------------------------------------------------------------------------- 
    927 %        \zstar- or \sstar-coordinate 
    928 % ------------------------------------------------------------------------------------------------------------- 
    929 \subsection{$Z^*$- or $S^*$-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.}) } 
     946%        z*- or s*-coordinate 
     947% ------------------------------------------------------------------------------------------------------------- 
     948\subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 
    930949\label{subsec:DOM_zgr_star} 
    931950 
    932 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.  
     951This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 
    933952 
    934953%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances 
     
    940959\label{subsec:DOM_msk} 
    941960 
    942 Whatever the vertical coordinate used, 
    943 the model offers the possibility of representing the bottom topography with steps that 
    944 follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 
    945 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, 
    946 which gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 
    947 mbathy is computed from the meter bathymetry using the definiton of gdept as 
    948 the number of $t$-points which gdept $\leq$ bathy. 
     961Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 
     962steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 
     963The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 
     964gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 
     965mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 
     966gdept $\leq$ bathy. 
    949967 
    950968Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after 
     
    954972As for the representation of bathymetry, a 2D integer array, misfdep, is created. 
    955973misfdep defines the level of the first wet $t$-point. 
    956 All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. 
    957 By default, misfdep(:,:)=1 and no cells are masked. 
     974All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 
     975By default, $misfdep(:,:) = 1$ and no cells are masked. 
    958976 
    959977In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into  
    960978the cavities are performed in the \textit{zgr\_isf} routine. 
    961 The compatibility between ice shelf draft and bathymetry is checked.  
    962 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked).  
     979The compatibility between ice shelf draft and bathymetry is checked. 
     980All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 
    963981If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 
    964982the bathymetry or the ice shelf draft is dug to fit this constrain. 
    965 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.\\  
     983If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 
    966984 
    967985From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 
    968 \begin{align*} 
    969   tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\ 
    970     \; 1&   \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$  }    \\ 
    971     \; 0&   \text{ if $k > mbathy(i,j)$  }    \end{cases}     \\ 
    972   umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
    973   vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\ 
    974   fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
    975                &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
    976   wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1) 
    977 \end{align*} 
     986\begin{alignat*}{2} 
     987  tmask(i,j,k) &= &  & 
     988    \begin{cases} 
     989                  0 &\text{if $                  k  <    misfdep(i,j)$} \\ 
     990                  1 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\ 
     991                  0 &\text{if $                  k  >     mbathy(i,j)$} 
     992    \end{cases} 
     993  \\ 
     994  umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     995  vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
     996  fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     997               &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     998  wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
     999  \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 
     1000\end{alignat*} 
    9781001 
    9791002Note that, without ice shelves cavities, 
    9801003masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 
    9811004Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface)  
    982 exactly in the same way as for the bottom boundary.  
     1005exactly in the same way as for the bottom boundary. 
    9831006 
    9841007The specification of closed lateral boundaries requires that at least 
    9851008the first and last rows and columns of the \textit{mbathy} array are set to zero. 
    986 In the particular case of an east-west cyclical boundary condition, 
    987 \textit{mbathy} has its last column equal to the second one and its first column equal to the last but one  
    988 (and so too the mask arrays) (see \autoref{fig:LBC_jperio}). 
    989  
     1009In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 
     1010the second one and its first column equal to the last but one (and so too the mask arrays) 
     1011(see \autoref{fig:LBC_jperio}). 
    9901012 
    9911013% ================================================================ 
     
    10001022 
    10011023Options are defined in \ngn{namtsd}. 
    1002 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 
     1024By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 
     1025salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 
    10031026\begin{description} 
    1004 \item[\np{ln\_tsd\_init}\forcode{ = .true.}] 
     1027\item[\np{ln\_tsd\_init}~\forcode{= .true.}] 
    10051028  use a T and S input files that can be given on the model grid itself or on their native input data grid. 
    10061029  In the latter case, 
     
    10091032  The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 
    10101033  The computation is done in the \mdl{dtatsd} module. 
    1011 \item[\np{ln\_tsd\_init}\forcode{ = .false.}] 
    1012   use constant salinity value of 35.5 psu and an analytical profile of temperature (typical of the tropical ocean), 
    1013   see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 
     1034\item[\np{ln\_tsd\_init}~\forcode{= .false.}] 
     1035  use constant salinity value of $35.5~psu$ and an analytical profile of temperature 
     1036  (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 
    10141037\end{description} 
    10151038 
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_TRA.tex

    r10442 r10502  
    1515%STEVEN :  is the use of the word "positive" to describe a scheme enough, or should it be "positive definite"? I added a comment to this effect on some instances of this below 
    1616 
    17 %\newpage 
    18  
    19 Using the representation described in \autoref{chap:DOM}, 
    20 several semi-discrete space forms of the tracer equations are available depending on 
    21 the vertical coordinate used and on the physics used. 
     17Using the representation described in \autoref{chap:DOM}, several semi -discrete space forms of 
     18the tracer equations are available depending on the vertical coordinate used and on the physics used. 
    2219In all the equations presented here, the masking has been omitted for simplicity. 
    23 One must be aware that all the quantities are masked fields and 
    24 that each time a mean or difference operator is used, 
    25 the resulting field is multiplied by a mask. 
     20One must be aware that all the quantities are masked fields and that each time a mean or 
     21difference operator is used, the resulting field is multiplied by a mask. 
    2622 
    2723The two active tracers are potential temperature and salinity. 
    2824Their prognostic equations can be summarized as follows: 
    2925\[ 
    30   \text{NXT} = \text{ADV}+\text{LDF}+\text{ZDF}+\text{SBC} 
    31   \ (+\text{QSR})\ (+\text{BBC})\ (+\text{BBL})\ (+\text{DMP}) 
     26  \text{NXT} =     \text{ADV} + \text{LDF} + \text{ZDF} + \text{SBC} 
     27               + \{\text{QSR},  \text{BBC},  \text{BBL},  \text{DMP}\} 
    3228\] 
    3329 
     
    3935The terms QSR, BBC, BBL and DMP are optional. 
    4036The external forcings and parameterisations require complex inputs and complex calculations 
    41 (\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, LDF and ZDF modules and  
    42 described in \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 
    43 Note that \mdl{tranpc}, the non-penetrative convection module, although located in the NEMO/OPA/TRA directory as 
    44 it directly modifies the tracer fields, is described with the model vertical physics (ZDF) together with 
     37(\eg bulk formulae, estimation of mixing coefficients) that are carried out in the SBC, 
     38LDF and ZDF modules and described in \autoref{chap:SBC}, \autoref{chap:LDF} and 
     39\autoref{chap:ZDF}, respectively. 
     40Note that \mdl{tranpc}, the non-penetrative convection module, although located in 
     41the \path{./src/OCE/TRA} directory as it directly modifies the tracer fields, 
     42is described with the model vertical physics (ZDF) together with 
    4543other available parameterization of convection. 
    4644 
     
    5048 
    5149The different options available to the user are managed by namelist logicals or CPP keys. 
    52 For each equation term  \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 
     50For each equation term \textit{TTT}, the namelist logicals are \textit{ln\_traTTT\_xxx}, 
    5351where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 
    5452The CPP key (when it exists) is \key{traTTT}. 
    5553The equivalent code can be found in the \textit{traTTT} or \textit{traTTT\_xxx} module, 
    56 in the NEMO/OPA/TRA directory. 
     54in the \path{./src/OCE/TRA} directory. 
    5755 
    5856The user has the option of extracting each tendency term on the RHS of the tracer equation for output 
    59 (\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}\forcode{ = .true.}), as described in \autoref{chap:DIA}. 
     57(\np{ln\_tra\_trd} or \np{ln\_tra\_mxl}~\forcode{= .true.}), as described in \autoref{chap:DIA}. 
    6058 
    6159% ================================================================ 
     
    7573\begin{equation} 
    7674  \label{eq:tra_adv} 
    77   ADV_\tau =-\frac{1}{b_t} \left( 
    78     \;\delta_i \left[ e_{2u}\,e_{3u} \;  u\; \tau_u  \right] 
    79     +\delta_j \left[ e_{1v}\,e_{3v}  \;  v\; \tau_v  \right] \; \right) 
    80   -\frac{1}{e_{3t}} \;\delta_k \left[ w\; \tau_w \right] 
     75  ADV_\tau = - \frac{1}{b_t} \Big(   \delta_i [ e_{2u} \, e_{3u} \; u \; \tau_u] 
     76                                   + \delta_j [ e_{1v} \, e_{3v} \; v \; \tau_v] \Big) 
     77             - \frac{1}{e_{3t}} \delta_k [w \; \tau_w] 
    8178\end{equation} 
    82 where $\tau$ is either T or S, and $b_t= e_{1t}\,e_{2t}\,e_{3t}$ is the volume of $T$-cells. 
     79where $\tau$ is either T or S, and $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells. 
    8380The flux form in \autoref{eq:tra_adv} implicitly requires the use of the continuity equation. 
    84 Indeed, it is obtained by using the following equality: 
    85 $\nabla \cdot \left( \vect{U}\,T \right)=\vect{U} \cdot \nabla T$ which 
    86 results from the use of the continuity equation,  $\partial _t e_3 + e_3\;\nabla \cdot \vect{U}=0$ 
    87 (which reduces to $\nabla \cdot \vect{U}=0$ in linear free surface, \ie \np{ln\_linssh}\forcode{ = .true.}). 
     81Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which 
     82results from the use of the continuity equation, $\partial_t e_3 + e_3 \; \nabla \cdot \vect U = 0$ 
     83(which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie \np{ln\_linssh}~\forcode{= .true.}). 
    8884Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that 
    8985it is consistent with the continuity equation in order to enforce the conservation properties of 
     
    9490\begin{figure}[!t] 
    9591  \begin{center} 
    96     \includegraphics[width=0.9\textwidth]{Fig_adv_scheme} 
     92    \includegraphics[]{Fig_adv_scheme} 
    9793    \caption{ 
    9894      \protect\label{fig:adv_scheme} 
     
    120116since the normal velocity is zero there. 
    121117At the sea surface the boundary condition depends on the type of sea surface chosen: 
     118 
    122119\begin{description} 
    123120\item[linear free surface:] 
    124   (\np{ln\_linssh}\forcode{ = .true.}) 
     121  (\np{ln\_linssh}~\forcode{= .true.}) 
    125122  the first level thickness is constant in time: 
    126   the vertical boundary condition is applied at the fixed surface $z=0$ rather than on the moving surface $z=\eta$. 
     123  the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on 
     124  the moving surface $z = \eta$. 
    127125  There is a non-zero advective flux which is set for all advection schemes as 
    128   $\left. {\tau_w } \right|_{k=1/2} =T_{k=1} $, 
    129   \ie the product of surface velocity (at $z=0$) by the first level tracer value. 
     126  $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie the product of surface velocity (at $z = 0$) by 
     127  the first level tracer value. 
    130128\item[non-linear free surface:] 
    131   (\np{ln\_linssh}\forcode{ = .false.}) 
     129  (\np{ln\_linssh}~\forcode{= .false.}) 
    132130  convergence/divergence in the first ocean level moves the free surface up/down. 
    133131  There is no tracer advection through it so that the advective fluxes through the surface are also zero. 
    134132\end{description} 
     133 
    135134In all cases, this boundary condition retains local conservation of tracer. 
    136135Global conservation is obtained in non-linear free surface case, but \textit{not} in the linear free surface case. 
    137136Nevertheless, in the latter case, it is achieved to a good approximation since 
    138137the non-conservative term is the product of the time derivative of the tracer and the free surface height, 
    139 two quantities that are not correlated \citep{Roullet_Madec_JGR00,Griffies_al_MWR01,Campin2004}. 
    140  
    141 The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) 
    142 is the centred (\textit{now}) \textit{effective} ocean velocity, 
    143 \ie the \textit{eulerian} velocity (see \autoref{chap:DYN}) plus 
    144 the eddy induced velocity (\textit{eiv}) and/or 
    145 the mixed layer eddy induced velocity (\textit{eiv}) when 
    146 those parameterisations are used (see \autoref{chap:LDF}). 
     138two quantities that are not correlated \citep{Roullet_Madec_JGR00, Griffies_al_MWR01, Campin2004}. 
     139 
     140The velocity field that appears in (\autoref{eq:tra_adv} and \autoref{eq:tra_adv_zco}) is 
     141the centred (\textit{now}) \textit{effective} ocean velocity, \ie the \textit{eulerian} velocity 
     142(see \autoref{chap:DYN}) plus the eddy induced velocity (\textit{eiv}) and/or 
     143the mixed layer eddy induced velocity (\textit{eiv}) when those parameterisations are used 
     144(see \autoref{chap:LDF}). 
    147145 
    148146Several tracer advection scheme are proposed, namely a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 
    149 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), 
    150 a Monotone Upstream Scheme for Conservative Laws scheme (MUSCL), 
    151 a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 
    152 and a Quadratic Upstream Interpolation for Convective Kinematics with 
    153 Estimated Streaming Terms scheme (QUICKEST). 
    154 The choice is made in the \textit{\ngn{namtra\_adv}} namelist, 
    155 by setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}. 
    156 The corresponding code can be found in the \mdl{traadv\_xxx} module, 
    157 where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 
    158 By default (\ie in the reference namelist, \ngn{namelist\_ref}), all the logicals are set to \forcode{.false.}. 
    159 If the user does not select an advection scheme in the configuration namelist (\ngn{namelist\_cfg}), 
     147a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), a Monotone Upstream Scheme for 
     148Conservative Laws scheme (MUSCL), a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 
     149and a Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms scheme (QUICKEST). 
     150The choice is made in the \ngn{namtra\_adv} namelist, by setting to \forcode{.true.} one of 
     151the logicals \textit{ln\_traadv\_xxx}. 
     152The corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 
     153\textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 
     154By default (\ie in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}. 
     155If the user does not select an advection scheme in the configuration namelist (\textit{namelist\_cfg}), 
    160156the tracers will \textit{not} be advected! 
    161157 
     
    163159The choosing an advection scheme is a complex matter which depends on the model physics, model resolution, 
    164160type of tracer, as well as the issue of numerical cost. In particular, we note that 
    165 (1) CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 
    166 they do not necessarily need additional diffusion; 
    167 (2) CEN and UBS are not \textit{positive} schemes 
    168 \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 
    169 implying that false extrema are permitted. 
    170 Their use is not recommended on passive tracers; 
    171 (3) It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 
    172 Indeed, if a source or sink of a passive tracer depends on an active one, 
    173 the difference of treatment of active and passive tracers can create very nice-looking frontal structures that 
    174 are pure numerical artefacts. 
     161 
     162\begin{enumerate} 
     163\item 
     164  CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 
     165  they do not necessarily need additional diffusion; 
     166\item 
     167  CEN and UBS are not \textit{positive} schemes 
     168  \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 
     169  implying that false extrema are permitted. 
     170  Their use is not recommended on passive tracers; 
     171\item 
     172  It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 
     173\end{enumerate} 
     174 
     175Indeed, if a source or sink of a passive tracer depends on an active one, the difference of treatment of active and 
     176passive tracers can create very nice-looking frontal structures that are pure numerical artefacts. 
    175177Nevertheless, most of our users set a different treatment on passive and active tracers, 
    176178that's the reason why this possibility is offered. 
    177 We strongly suggest them to perform a sensitivity experiment using a same treatment to 
    178 assess the robustness of their results. 
     179We strongly suggest them to perform a sensitivity experiment using a same treatment to assess the robustness of 
     180their results. 
    179181 
    180182% ------------------------------------------------------------------------------------------------------------- 
    181183%        2nd and 4th order centred schemes 
    182184% ------------------------------------------------------------------------------------------------------------- 
    183 \subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}\forcode{ = .true.})} 
     185\subsection{CEN: Centred scheme (\protect\np{ln\_traadv\_cen}~\forcode{= .true.})} 
    184186\label{subsec:TRA_adv_cen} 
    185187 
    186188%        2nd order centred scheme   
    187189 
    188 The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}\forcode{ = .true.}. 
     190The centred advection scheme (CEN) is used when \np{ln\_traadv\_cen}~\forcode{= .true.}. 
    189191Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 
    190192setting \np{nn\_cen\_h} and \np{nn\_cen\_v} to $2$ or $4$. 
     
    196198\begin{equation} 
    197199  \label{eq:tra_adv_cen2} 
    198   \tau_u^{cen2} =\overline T ^{i+1/2} 
     200  \tau_u^{cen2} = \overline T ^{i + 1/2} 
    199201\end{equation} 
    200202 
    201 CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2)$ but dispersive 
     203CEN2 is non diffusive (\ie it conserves the tracer variance, $\tau^2$) but dispersive 
    202204(\ie it may create false extrema). 
    203205It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 
    204206produce a sensible solution. 
    205207The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 
    206 so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value.  
     208so $T$ in (\autoref{eq:tra_adv_cen2}) is the \textit{now} tracer value. 
    207209 
    208210Note that using the CEN2, the overall tracer advection is of second order accuracy since 
     
    216218\begin{equation} 
    217219  \label{eq:tra_adv_cen4} 
    218   \tau_u^{cen4} =\overline{   T - \frac{1}{6}\,\delta_i \left[ \delta_{i+1/2}[T] \,\right]   }^{\,i+1/2} 
     220  \tau_u^{cen4} = \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \, \Big]}^{\,i + 1/2} 
    219221\end{equation} 
    220 In the vertical direction (\np{nn\_cen\_v}\forcode{ = 4}), 
     222In the vertical direction (\np{nn\_cen\_v}~\forcode{= 4}), 
    221223a $4^{th}$ COMPACT interpolation has been prefered \citep{Demange_PhD2014}. 
    222224In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion, 
    223 spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}. 
    224   
     225spectral characteristics similar to schemes of higher order \citep{Lele_JCP1992}.  
    225226 
    226227Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but 
     
    249250%        FCT scheme   
    250251% ------------------------------------------------------------------------------------------------------------- 
    251 \subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}\forcode{ = .true.})} 
     252\subsection{FCT: Flux Corrected Transport scheme (\protect\np{ln\_traadv\_fct}~\forcode{= .true.})} 
    252253\label{subsec:TRA_adv_tvd} 
    253254 
    254 The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}\forcode{ = .true.}. 
     255The Flux Corrected Transport schemes (FCT) is used when \np{ln\_traadv\_fct}~\forcode{= .true.}. 
    255256Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 
    256257setting \np{nn\_fct\_h} and \np{nn\_fct\_v} to $2$ or $4$. 
     
    263264  \label{eq:tra_adv_fct} 
    264265  \begin{split} 
    265     \tau_u^{ups}&= 
     266    \tau_u^{ups} &= 
    266267    \begin{cases} 
    267       T_{i+1}  & \text{if $\ u_{i+1/2} <     0$} \hfill \\ 
    268       T_i         & \text{if $\ u_{i+1/2} \geq 0$} \hfill \\ 
     268                     T_{i + 1} & \text{if~} u_{i + 1/2} <    0 \\ 
     269                     T_i       & \text{if~} u_{i + 1/2} \geq 0 \\ 
    269270    \end{cases} 
    270     \\ \\ 
    271     \tau_u^{fct}&=\tau_u^{ups} +c_u \;\left( {\tau_u^{cen} -\tau_u^{ups} } \right) 
     271    \\ 
     272    \tau_u^{fct} &= \tau_u^{ups} + c_u \, \big( \tau_u^{cen} - \tau_u^{ups} \big) 
    272273  \end{split} 
    273274\end{equation} 
     
    278279The one chosen in \NEMO is described in \citet{Zalesak_JCP79}. 
    279280$c_u$ only departs from $1$ when the advective term produces a local extremum in the tracer field. 
    280 The resulting scheme is quite expensive but \emph{positive}. 
     281The resulting scheme is quite expensive but \textit{positive}. 
    281282It can be used on both active and passive tracers. 
    282283A comparison of FCT-2 with MUSCL and a MPDATA scheme can be found in \citet{Levy_al_GRL01}. 
     
    294295$\tau_u^{ups}$ is evaluated using the \textit{before} tracer. 
    295296In other words, the advective part of the scheme is time stepped with a leap-frog scheme 
    296 while a forward scheme is used for the diffusive part.  
     297while a forward scheme is used for the diffusive part. 
    297298 
    298299% ------------------------------------------------------------------------------------------------------------- 
    299300%        MUSCL scheme   
    300301% ------------------------------------------------------------------------------------------------------------- 
    301 \subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}\forcode{ = .true.})} 
     302\subsection{MUSCL: Monotone Upstream Scheme for Conservative Laws (\protect\np{ln\_traadv\_mus}~\forcode{= .true.})} 
    302303\label{subsec:TRA_adv_mus} 
    303304 
    304 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}\forcode{ = .true.}. 
     305The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np{ln\_traadv\_mus}~\forcode{= .true.}. 
    305306MUSCL implementation can be found in the \mdl{traadv\_mus} module. 
    306307 
     
    309310two $T$-points (\autoref{fig:adv_scheme}). 
    310311For example, in the $i$-direction : 
    311 \[ 
     312\begin{equation} 
    312313  % \label{eq:tra_adv_mus} 
    313   \tau_u^{mus} = \left\{ 
    314     \begin{aligned} 
    315       &\tau_i  &+ \frac{1}{2} \;\left( 1-\frac{u_{i+1/2} \;\rdt}{e_{1u}} \right) 
    316       &\ \widetilde{\partial _i \tau}  & \quad \text{if }\;u_{i+1/2} \geqslant 0      \\ 
    317       &\tau_{i+1/2} &+\frac{1}{2}\;\left( 1+\frac{u_{i+1/2} \;\rdt}{e_{1u} } \right) 
    318       &\ \widetilde{\partial_{i+1/2} \tau } & \text{if }\;u_{i+1/2} <0 
    319     \end{aligned} 
    320   \right. 
    321 \] 
    322 where $\widetilde{\partial _i \tau}$ is the slope of the tracer on which a limitation is imposed to 
     314  \tau_u^{mus} = \lt\{ 
     315  \begin{split} 
     316                       \tau_i         &+ \frac{1}{2} \lt( 1 - \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
     317                       \widetilde{\partial_i         \tau} & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
     318                       \tau_{i + 1/2} &+ \frac{1}{2} \lt( 1 + \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
     319                       \widetilde{\partial_{i + 1/2} \tau} & \text{if~} u_{i + 1/2} <         0 
     320  \end{split} 
     321                                                                                                      \rt. 
     322\end{equation} 
     323where $\widetilde{\partial_i \tau}$ is the slope of the tracer on which a limitation is imposed to 
    323324ensure the \textit{positive} character of the scheme. 
    324325 
    325 The time stepping is performed using a forward scheme, 
    326 that is the \textit{before} tracer field is used to evaluate $\tau_u^{mus}$. 
     326The time stepping is performed using a forward scheme, that is the \textit{before} tracer field is used to 
     327evaluate $\tau_u^{mus}$. 
    327328 
    328329For an ocean grid point adjacent to land and where the ocean velocity is directed toward land, 
     
    330331This choice ensure the \textit{positive} character of the scheme. 
    331332In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes 
    332 (\np{ln\_mus\_ups}\forcode{ = .true.}). 
     333(\np{ln\_mus\_ups}~\forcode{= .true.}). 
    333334 
    334335% ------------------------------------------------------------------------------------------------------------- 
    335336%        UBS scheme   
    336337% ------------------------------------------------------------------------------------------------------------- 
    337 \subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}\forcode{ = .true.})} 
     338\subsection{UBS a.k.a. UP3: Upstream-Biased Scheme (\protect\np{ln\_traadv\_ubs}~\forcode{= .true.})} 
    338339\label{subsec:TRA_adv_ubs} 
    339340 
    340 The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}\forcode{ = .true.}. 
     341The Upstream-Biased Scheme (UBS) is used when \np{ln\_traadv\_ubs}~\forcode{= .true.}. 
    341342UBS implementation can be found in the \mdl{traadv\_mus} module. 
    342343 
     
    347348\begin{equation} 
    348349  \label{eq:tra_adv_ubs} 
    349   \tau_u^{ubs} =\overline T ^{i+1/2}-\;\frac{1}{6} \left\{ 
    350     \begin{aligned} 
    351       &\tau"_i          & \quad \text{if }\ u_{i+1/2} \geqslant 0      \\ 
    352       &\tau"_{i+1}   & \quad \text{if }\ u_{i+1/2}       <       0 
    353     \end{aligned} 
    354   \right. 
     350  \tau_u^{ubs} = \overline T ^{i + 1/2} - \frac{1}{6} 
     351    \begin{cases} 
     352                                                      \tau"_i       & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
     353                                                      \tau"_{i + 1} & \text{if~} u_{i + 1/2} <         0 
     354    \end{cases} 
     355  \quad 
     356  \text{where~} \tau"_i = \delta_i \lt[ \delta_{i + 1/2} [\tau] \rt] 
    355357\end{equation} 
    356 where $\tau "_i =\delta_i \left[ {\delta_{i+1/2} \left[ \tau \right]} \right]$. 
    357  
    358 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 
     358 
     359This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 
    359360\citep{Shchepetkin_McWilliams_OM05}. 
    360361The overall performance of the advection scheme is similar to that reported in \cite{Farrow1995}. 
    361362It is a relatively good compromise between accuracy and smoothness. 
    362 Nevertheless the scheme is not \emph{positive}, meaning that false extrema are permitted, 
     363Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted, 
    363364but the amplitude of such are significantly reduced over the centred second or fourth order method. 
    364365Therefore it is not recommended that it should be applied to a passive tracer that requires positivity. 
     
    368369\citep{Shchepetkin_McWilliams_OM05, Demange_PhD2014}. 
    369370Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme 
    370 (\np{nn\_cen\_v}\forcode{ = 2 or 4}). 
     371(\np{nn\_cen\_v}~\forcode{= 2 or 4}). 
    371372 
    372373For stability reasons (see \autoref{chap:STP}), the first term  in \autoref{eq:tra_adv_ubs} 
     
    382383 
    383384Note that it is straightforward to rewrite \autoref{eq:tra_adv_ubs} as follows: 
    384 \[ 
     385\begin{gather} 
    385386  \label{eq:traadv_ubs2} 
    386   \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} \left\{ 
    387     \begin{aligned} 
    388       & + \tau"_i       & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ 
    389       &  - \tau"_{i+1}     & \quad \text{if }\ u_{i+1/2}       <       0 
    390     \end{aligned} 
    391   \right. 
    392 \] 
    393 or equivalently  
    394 \[ 
     387  \tau_u^{ubs} = \tau_u^{cen4} + \frac{1}{12} 
     388    \begin{cases} 
     389      + \tau"_i       & \text{if} \ u_{i + 1/2} \geqslant 0 \\ 
     390      - \tau"_{i + 1} & \text{if} \ u_{i + 1/2} <         0 
     391    \end{cases} 
     392  \intertext{or equivalently} 
    395393  % \label{eq:traadv_ubs2b} 
    396   u_{i+1/2} \ \tau_u^{ubs} 
    397   =u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\delta_i\left[ \delta_{i+1/2}[T] \,\right] }^{\,i+1/2} 
    398   - \frac{1}{2} |u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] 
    399 \] 
     394  u_{i + 1/2} \ \tau_u^{ubs} = u_{i + 1/2} \, \overline{T - \frac{1}{6} \, \delta_i \Big[ \delta_{i + 1/2}[T] \Big]}^{\,i + 1/2} 
     395                             - \frac{1}{2} |u|_{i + 1/2} \, \frac{1}{6} \, \delta_{i + 1/2} [\tau"_i] \nonumber 
     396\end{gather} 
    400397 
    401398\autoref{eq:traadv_ubs2} has several advantages. 
     
    403400an upstream-biased diffusion term is added. 
    404401Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to 
    405 be evaluated at the \emph{now} time step using \autoref{eq:tra_adv_ubs}. 
     402be evaluated at the \textit{now} time step using \autoref{eq:tra_adv_ubs}. 
    406403Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which 
    407 is simply proportional to the velocity: 
    408 $A_u^{lm}= \frac{1}{12}\,{e_{1u}}^3\,|u|$. 
     404is simply proportional to the velocity: $A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$. 
    409405Note the current version of NEMO uses the computationally more efficient formulation \autoref{eq:tra_adv_ubs}. 
    410406 
     
    412408%        QCK scheme   
    413409% ------------------------------------------------------------------------------------------------------------- 
    414 \subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}\forcode{ = .true.})} 
     410\subsection{QCK: QuiCKest scheme (\protect\np{ln\_traadv\_qck}~\forcode{= .true.})} 
    415411\label{subsec:TRA_adv_qck} 
    416412 
    417413The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme 
    418 proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}\forcode{ = .true.}. 
     414proposed by \citet{Leonard1979} is used when \np{ln\_traadv\_qck}~\forcode{= .true.}. 
    419415QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 
    420416 
     
    422418\citep{Leonard1991}. 
    423419It has been implemented in NEMO by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 
    424 The resulting scheme is quite expensive but \emph{positive}. 
     420The resulting scheme is quite expensive but \textit{positive}. 
    425421It can be used on both active and passive tracers. 
    426422However, the intrinsic diffusion of QCK makes its use risky in the vertical direction where 
     
    431427 
    432428%%%gmcomment   :  Cross term are missing in the current implementation.... 
    433  
    434429 
    435430% ================================================================ 
     
    458453except for the pure vertical component that appears when a rotation tensor is used. 
    459454This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:STP}). 
    460 When \np{ln\_traldf\_msc}\forcode{ = .true.}, a Method of Stabilizing Correction is used in which 
     455When \np{ln\_traldf\_msc}~\forcode{= .true.}, a Method of Stabilizing Correction is used in which 
    461456the pure vertical component is split into an explicit and an implicit part \citep{Lemarie_OM2012}. 
    462457 
     
    464459%        Type of operator 
    465460% ------------------------------------------------------------------------------------------------------------- 
    466 \subsection[Type of operator (\protect\np{ln\_traldf}\{\_NONE,\_lap,\_blp\}\})] 
    467               {Type of operator (\protect\np{ln\_traldf\_NONE}, \protect\np{ln\_traldf\_lap}, or \protect\np{ln\_traldf\_blp}) }  
     461\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}) }  
    468462\label{subsec:TRA_ldf_op} 
    469463 
    470464Three operator options are proposed and, one and only one of them must be selected: 
     465 
    471466\begin{description} 
    472 \item[\np{ln\_traldf\_NONE}\forcode{ = .true.}:] 
     467\item[\np{ln\_traldf\_NONE}~\forcode{= .true.}:] 
    473468  no operator selected, the lateral diffusive tendency will not be applied to the tracer equation. 
    474469  This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example). 
    475 \item[\np{ln\_traldf\_lap}\forcode{ = .true.}:] 
     470\item[\np{ln\_traldf\_lap}~\forcode{= .true.}:] 
    476471  a laplacian operator is selected. 
    477   This harmonic operator takes the following expression:  $\mathpzc{L}(T)=\nabla \cdot A_{ht}\;\nabla T $, 
     472  This harmonic operator takes the following expression:  $\mathpzc{L}(T) = \nabla \cdot A_{ht} \; \nabla T $, 
    478473  where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 
    479474  and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 
    480 \item[\np{ln\_traldf\_blp}\forcode{ = .true.}]: 
     475\item[\np{ln\_traldf\_blp}~\forcode{= .true.}]: 
    481476  a bilaplacian operator is selected. 
    482477  This biharmonic operator takes the following expression: 
    483   $\mathpzc{B}=- \mathpzc{L}\left(\mathpzc{L}(T) \right) = -\nabla \cdot b\nabla \left( {\nabla \cdot b\nabla T} \right)$ 
     478  $\mathpzc{B} = - \mathpzc{L}(\mathpzc{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$ 
    484479  where the gradient operats along the selected direction, 
    485   and $b^2=B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 
     480  and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 
    486481  In the code, the bilaplacian operator is obtained by calling the laplacian twice. 
    487482\end{description} 
     
    495490whereas the laplacian damping time scales only like $\lambda^{-2}$. 
    496491 
    497  
    498492% ------------------------------------------------------------------------------------------------------------- 
    499493%        Direction of action 
    500494% ------------------------------------------------------------------------------------------------------------- 
    501 \subsection[Action direction (\protect\np{ln\_traldf}\{\_lev,\_hor,\_iso,\_triad\})] 
    502               {Direction of action (\protect\np{ln\_traldf\_lev}, \protect\np{ln\_traldf\_hor}, \protect\np{ln\_traldf\_iso}, or \protect\np{ln\_traldf\_triad}) }  
     495\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}) }  
    503496\label{subsec:TRA_ldf_dir} 
    504497 
    505498The choice of a direction of action determines the form of operator used. 
    506499The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when 
    507 iso-level option is used (\np{ln\_traldf\_lev}\forcode{ = .true.}) or 
    508 when a horizontal (\ie geopotential) operator is demanded in \zstar-coordinate 
     500iso-level option is used (\np{ln\_traldf\_lev}~\forcode{= .true.}) or 
     501when a horizontal (\ie geopotential) operator is demanded in \textit{z}-coordinate 
    509502(\np{ln\_traldf\_hor} and \np{ln\_zco} equal \forcode{.true.}). 
    510503The associated code can be found in the \mdl{traldf\_lap\_blp} module. 
     
    521514 
    522515The resulting discret form of the three operators (one iso-level and two rotated one) is given in 
    523 the next two sub-sections.  
    524  
     516the next two sub-sections. 
    525517 
    526518% ------------------------------------------------------------------------------------------------------------- 
    527519%       iso-level operator 
    528520% ------------------------------------------------------------------------------------------------------------- 
    529 \subsection{Iso-level (bi-)laplacian operator ( \protect\np{ln\_traldf\_iso}) } 
     521\subsection{Iso-level (bi -)laplacian operator ( \protect\np{ln\_traldf\_iso}) } 
    530522\label{subsec:TRA_ldf_lev} 
    531523 
     
    533525\begin{equation} 
    534526  \label{eq:tra_ldf_lap} 
    535   D_t^{lT} =\frac{1}{b_t} \left( \; 
    536     \delta_{i}\left[ A_u^{lT} \; \frac{e_{2u}\,e_{3u}}{e_{1u}} \;\delta_{i+1/2} [T] \right] 
    537     + \delta_{j}\left[ A_v^{lT} \;  \frac{e_{1v}\,e_{3v}}{e_{2v}} \;\delta_{j+1/2} [T] \right]  \;\right) 
     527  D_t^{lT} = \frac{1}{b_t} \Bigg(   \delta_{i} \lt[ A_u^{lT} \; \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [T] \rt] 
     528                                  + \delta_{j} \lt[ A_v^{lT} \; \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [T] \rt] \Bigg) 
    538529\end{equation} 
    539 where  $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells and 
     530where $b_t = e_{1t} \, e_{2t} \, e_{3t}$  is the volume of $T$-cells and 
    540531where zero diffusive fluxes is assumed across solid boundaries, 
    541532first (and third in bilaplacian case) horizontal tracer derivative are masked. 
    542533It is implemented in the \rou{traldf\_lap} subroutine found in the \mdl{traldf\_lap} module. 
    543534The module also contains \rou{traldf\_blp}, the subroutine calling twice \rou{traldf\_lap} in order to 
    544 compute the iso-level bilaplacian operator.  
    545  
    546 It is a \emph{horizontal} operator (\ie acting along geopotential surfaces) in 
     535compute the iso-level bilaplacian operator. 
     536 
     537It is a \textit{horizontal} operator (\ie acting along geopotential surfaces) in 
    547538the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 
    548 It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}\forcode{ = .true.}, 
    549 we have \np{ln\_traldf\_lev}\forcode{ = .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}\forcode{ = .true.}. 
     539It is thus used when, in addition to \np{ln\_traldf\_lap} or \np{ln\_traldf\_blp}~\forcode{= .true.}, 
     540we have \np{ln\_traldf\_lev}~\forcode{= .true.} or \np{ln\_traldf\_hor}~=~\np{ln\_zco}~\forcode{= .true.}. 
    550541In both cases, it significantly contributes to diapycnal mixing. 
    551542It is therefore never recommended, even when using it in the bilaplacian case. 
    552543 
    553 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), 
     544Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}), 
    554545tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom. 
    555546In this case, horizontal derivatives in (\autoref{eq:tra_ldf_lap}) at the bottom level require a specific treatment. 
    556547They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}. 
    557548 
    558  
    559549% ------------------------------------------------------------------------------------------------------------- 
    560550%         Rotated laplacian operator 
    561551% ------------------------------------------------------------------------------------------------------------- 
    562 \subsection{Standard and triad (bi-)laplacian operator} 
     552\subsection{Standard and triad (bi -)laplacian operator} 
    563553\label{subsec:TRA_ldf_iso_triad} 
    564554 
    565 %&&    Standard rotated (bi-)laplacian operator 
     555%&&    Standard rotated (bi -)laplacian operator 
    566556%&& ---------------------------------------------- 
    567 \subsubsection{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})} 
     557\subsubsection{Standard rotated (bi -)laplacian operator (\protect\mdl{traldf\_iso})} 
    568558\label{subsec:TRA_ldf_iso} 
    569  
    570559The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:PE_zdf}) 
    571 takes the following semi-discrete space form in $z$- and $s$-coordinates: 
     560takes the following semi -discrete space form in $z$- and $s$-coordinates: 
    572561\begin{equation} 
    573562  \label{eq:tra_ldf_iso} 
    574563  \begin{split} 
    575     D_T^{lT} = \frac{1}{b_t}   & \left\{   \,\;\delta_i \left[   A_u^{lT}  \left( 
    576           \frac{e_{2u}\,e_{3u}}{e_{1u}} \,\delta_{i+1/2}[T] 
    577           - e_{2u}\;r_{1u} \,\overline{\overline{ \delta_{k+1/2}[T] }}^{\,i+1/2,k} 
    578         \right)   \right]   \right.    \\ 
    579     &             +\delta_j \left[ A_v^{lT} \left( 
    580         \frac{e_{1v}\,e_{3v}}{e_{2v}}  \,\delta_{j+1/2} [T] 
    581         - e_{1v}\,r_{2v} \,\overline{\overline{ \delta_{k+1/2} [T] }}^{\,j+1/2,k} 
    582       \right)   \right]                 \\ 
    583     & +\delta_k \left[ A_w^{lT} \left( 
    584         -\;e_{2w}\,r_{1w} \,\overline{\overline{ \delta_{i+1/2} [T] }}^{\,i,k+1/2} 
    585       \right.   \right.                 \\ 
    586     & \qquad \qquad \quad 
    587     - e_{1w}\,r_{2w} \,\overline{\overline{ \delta_{j+1/2} [T] }}^{\,j,k+1/2}     \\ 
    588     & \left. {\left. {   \qquad \qquad \ \ \ \left. { 
    589                 +\;\frac{e_{1w}\,e_{2w}}{e_{3w}} \,\left( r_{1w}^2 + r_{2w}^2 \right) 
    590                 \,\delta_{k+1/2} [T] } \right) } \right] \quad } \right\} 
     564    D_T^{lT} = \frac{1}{b_t} \Bigg[ \quad &\delta_i A_u^{lT} \lt( \frac{e_{2u} e_{3u}}{e_{1u}}                      \, \delta_{i + 1/2} [T] 
     565                                                                  - e_{2u} r_{1u} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,i + 1/2,k} \rt) \Bigg. \\ 
     566                                    +     &\delta_j A_v^{lT} \lt( \frac{e_{1v} e_{3v}}{e_{2v}}                       \, \delta_{j + 1/2} [T] 
     567                                                                  - e_{1v} r_{2v} \, \overline{\overline{\delta_{k + 1/2} [T]}}^{\,j + 1/2,k} \rt)        \\ 
     568                                   +     &\delta_k A_w^{lT} \lt( \frac{e_{1w} e_{2w}}{e_{3w}} (r_{1w}^2 + r_{2w}^2) \, \delta_{k + 1/2} [T] \rt.           \\ 
     569                                   & \qquad \quad \Bigg. \lt.     - e_{2w} r_{1w} \, \overline{\overline{\delta_{i + 1/2} [T]}}^{\,i,k + 1/2} 
     570                                                                  - e_{1w} r_{2w} \, \overline{\overline{\delta_{j + 1/2} [T]}}^{\,j,k + 1/2} \rt) \Bigg] 
    591571  \end{split} 
    592572\end{equation} 
    593 where $b_t$=$e_{1t}\,e_{2t}\,e_{3t}$  is the volume of $T$-cells, 
     573where $b_t = e_{1t} \, e_{2t} \, e_{3t}$  is the volume of $T$-cells, 
    594574$r_1$ and $r_2$ are the slopes between the surface of computation ($z$- or $s$-surfaces) and 
    595575the surface along which the diffusion operator acts (\ie horizontal or iso-neutral surfaces). 
    596 It is thus used when, in addition to \np{ln\_traldf\_lap}\forcode{ = .true.}, 
    597 we have \np{ln\_traldf\_iso}\forcode{ = .true.}, 
    598 or both \np{ln\_traldf\_hor}\forcode{ = .true.} and \np{ln\_zco}\forcode{ = .true.}. 
     576It is thus used when, in addition to \np{ln\_traldf\_lap}~\forcode{= .true.}, 
     577we have \np{ln\_traldf\_iso}~\forcode{= .true.}, 
     578or both \np{ln\_traldf\_hor}~\forcode{= .true.} and \np{ln\_zco}~\forcode{= .true.}. 
    599579The way these slopes are evaluated is given in \autoref{sec:LDF_slp}. 
    600580At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using 
    601 the mask technique (see \autoref{sec:LBC_coast}).  
     581the mask technique (see \autoref{sec:LBC_coast}). 
    602582 
    603583The operator in \autoref{eq:tra_ldf_iso} involves both lateral and vertical derivatives. 
     
    606586For computer efficiency reasons, this term is not computed in the \mdl{traldf\_iso} module, 
    607587but in the \mdl{trazdf} module where, if iso-neutral mixing is used, 
    608 the vertical mixing coefficient is simply increased by 
    609 $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$.  
     588the vertical mixing coefficient is simply increased by $\frac{e_{1w} e_{2w}}{e_{3w}}(r_{1w}^2 + r_{2w}^2)$. 
    610589 
    611590This formulation conserves the tracer but does not ensure the decrease of the tracer variance. 
    612591Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without 
    613 any additional background horizontal diffusion \citep{Guilyardi_al_CD01}.  
    614  
    615 Note that in the partial step $z$-coordinate (\np{ln\_zps}\forcode{ = .true.}), 
     592any additional background horizontal diffusion \citep{Guilyardi_al_CD01}. 
     593 
     594Note that in the partial step $z$-coordinate (\np{ln\_zps}~\forcode{= .true.}), 
    616595the horizontal derivatives at the bottom level in \autoref{eq:tra_ldf_iso} require a specific treatment. 
    617596They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}. 
    618597 
    619 %&&     Triad rotated (bi-)laplacian operator 
     598%&&     Triad rotated (bi -)laplacian operator 
    620599%&&  ------------------------------------------- 
    621 \subsubsection{Triad rotated (bi-)laplacian operator (\protect\np{ln\_traldf\_triad})} 
     600\subsubsection{Triad rotated (bi -)laplacian operator (\protect\np{ln\_traldf\_triad})} 
    622601\label{subsec:TRA_ldf_triad} 
    623602 
    624 If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}\forcode{ = .true.}; see \autoref{apdx:triad}) 
     603If the Griffies triad scheme is employed (\np{ln\_traldf\_triad}~\forcode{= .true.}; see \autoref{apdx:triad}) 
    625604 
    626605An alternative scheme developed by \cite{Griffies_al_JPO98} which ensures tracer variance decreases 
    627 is also available in \NEMO (\np{ln\_traldf\_grif}\forcode{ = .true.}). 
     606is also available in \NEMO (\np{ln\_traldf\_grif}~\forcode{= .true.}). 
    628607A complete description of the algorithm is given in \autoref{apdx:triad}. 
    629608 
     
    635614It requires an additional assumption on boundary conditions: 
    636615first and third derivative terms normal to the coast, 
    637 normal to the bottom and normal to the surface are set to zero.  
     616normal to the bottom and normal to the surface are set to zero. 
    638617 
    639618%&&    Option for the rotated operators 
     
    642621\label{subsec:TRA_ldf_options} 
    643622 
    644 \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 
    645  
    646 \np{rn\_slpmax} = slope limit (both operators) 
    647  
    648 \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 
    649  
    650 \np{rn\_sw\_triad} =1 switching triad; 
    651                    =0 all 4 triads used (triad only)  
    652  
    653 \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 
     623\begin{itemize} 
     624\item \np{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 
     625\item \np{rn\_slpmax} = slope limit (both operators) 
     626\item \np{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 
     627\item \np{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only)  
     628\item \np{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 
     629\end{itemize} 
    654630 
    655631% ================================================================ 
     
    666642The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, 
    667643and is based on a laplacian operator. 
    668 The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi-discrete space form: 
    669 \[ 
     644The vertical diffusion operator given by (\autoref{eq:PE_zdf}) takes the following semi -discrete space form: 
     645\begin{gather*} 
    670646  % \label{eq:tra_zdf} 
    671   \begin{split} 
    672     D^{vT}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vT}_w}{e_{3w}}  \delta_{k+1/2}[T] \;\right]    \\ 
    673     D^{vS}_T &= \frac{1}{e_{3t}} \; \delta_k \left[ \;\frac{A^{vS}_w}{e_{3w}}  \delta_{k+1/2}[S] \;\right] 
    674   \end{split} 
    675 \] 
     647    D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \\ 
     648    D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt] 
     649\end{gather*} 
    676650where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on temperature and salinity, 
    677651respectively. 
    678 Generally, $A_w^{vT}=A_w^{vS}$ except when double diffusive mixing is parameterised (\ie \key{zdfddm} is defined). 
     652Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised 
     653(\ie \key{zdfddm} is defined). 
    679654The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF). 
    680655Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by 
    681 $\frac{e_{1w}\,e_{2w} }{e_{3w} }\ \left( {r_{1w} ^2+r_{2w} ^2} \right)$ to account for 
    682 the vertical second derivative of \autoref{eq:tra_ldf_iso}. 
     656$\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to account for the vertical second derivative of 
     657\autoref{eq:tra_ldf_iso}. 
    683658 
    684659At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified. 
    685660At the surface they are prescribed from the surface forcing and added in a dedicated routine 
    686661(see \autoref{subsec:TRA_sbc}), whilst at the bottom they are set to zero for heat and salt unless 
    687 a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}).  
     662a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 
    688663 
    689664The large eddy coefficient found in the mixed layer together with high vertical resolution implies that 
    690 in the case of explicit time stepping (\np{ln\_zdfexp}\forcode{ = .true.}) 
     665in the case of explicit time stepping (\np{ln\_zdfexp}~\forcode{= .true.}) 
    691666there would be too restrictive a constraint on the time step. 
    692667Therefore, the default implicit time stepping is preferred for the vertical diffusion since 
    693668it overcomes the stability constraint. 
    694 A forward time differencing scheme (\np{ln\_zdfexp}\forcode{ = .true.}) using 
     669A forward time differencing scheme (\np{ln\_zdfexp}~\forcode{= .true.}) using 
    695670a time splitting technique (\np{nn\_zdfexp} $> 1$) is provided as an alternative. 
    696 Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics.  
     671Namelist variables \np{ln\_zdfexp} and \np{nn\_zdfexp} apply to both tracers and dynamics. 
    697672 
    698673% ================================================================ 
     
    712687This has been found to enhance readability of the code. 
    713688The two formulations are completely equivalent; 
    714 the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer.  
     689the forcing terms in trasbc are the surface fluxes divided by the thickness of the top model layer. 
    715690 
    716691Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 
     
    724699The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields (used on tracers): 
    725700 
    726 $\bullet$ $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 
    727 (\ie the difference between the total surface heat flux and the fraction of the short wave flux that  
    728 penetrates into the water column, see \autoref{subsec:TRA_qsr}) 
    729 plus the heat content associated with of the mass exchange with the atmosphere and lands. 
    730  
    731 $\bullet$ $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 
    732  
    733 $\bullet$ \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 
    734 possibly with the sea-ice and ice-shelves. 
    735  
    736 $\bullet$ \textit{rnf}, the mass flux associated with runoff 
    737 (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
    738  
    739 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 
    740 (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 
     701\begin{itemize} 
     702\item 
     703  $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 
     704  (\ie the difference between the total surface heat flux and the fraction of the short wave flux that 
     705  penetrates into the water column, see \autoref{subsec:TRA_qsr}) 
     706  plus the heat content associated with of the mass exchange with the atmosphere and lands. 
     707\item 
     708  $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 
     709\item 
     710  \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 
     711  possibly with the sea-ice and ice-shelves. 
     712\item 
     713  \textit{rnf}, the mass flux associated with runoff 
     714  (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
     715\item 
     716  \textit{fwfisf}, the mass flux associated with ice shelf melt, 
     717  (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 
     718\end{itemize} 
    741719 
    742720The surface boundary condition on temperature and salinity is applied as follows: 
    743721\begin{equation} 
    744722  \label{eq:tra_sbc} 
    745   \begin{aligned} 
    746     &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} }  &\overline{ Q_{ns}       }^t  & \\ 
    747     & F^S =\frac{ 1 }{\rho_o  \,      \left. e_{3t} \right|_{k=1} }  &\overline{ \textit{sfx} }^t   & \\ 
    748   \end{aligned} 
    749 \end{equation}  
    750 where $\overline{x }^t$ means that $x$ is averaged over two consecutive time steps ($t-\rdt/2$ and $t+\rdt/2$). 
     723  \begin{alignedat}{2} 
     724    F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{Q_{ns}      }^t \\ 
     725    F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{\textit{sfx}}^t 
     726  \end{alignedat} 
     727\end{equation} 
     728where $\overline x^t$ means that $x$ is averaged over two consecutive time steps 
     729($t - \rdt / 2$ and $t + \rdt / 2$). 
    751730Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:STP}). 
    752731 
    753 In the linear free surface case (\np{ln\_linssh}\forcode{ = .true.}), an additional term has to be added on 
     732In the linear free surface case (\np{ln\_linssh}~\forcode{= .true.}), an additional term has to be added on 
    754733both temperature and salinity. 
    755734On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$. 
     
    759738\begin{equation} 
    760739  \label{eq:tra_sbc_lin} 
    761   \begin{aligned} 
    762     &F^T = \frac{ 1 }{\rho_o \;C_p \,\left. e_{3t} \right|_{k=1} } 
    763     &\overline{ \left( Q_{ns} - \textit{emp}\;C_p\,\left. T \right|_{k=1} \right) }^t  & \\ 
    764     % 
    765     & F^S =\frac{ 1 }{\rho_o \,\left. e_{3t} \right|_{k=1} } 
    766     &\overline{ \left( \;\textit{sfx} - \textit{emp} \;\left. S \right|_{k=1}  \right) }^t   & \\ 
    767   \end{aligned} 
     740  \begin{alignedat}{2} 
     741    F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
     742          &\overline{(Q_{ns}       - C_p \, \textit{emp} \lt. T \rt|_{k = 1})}^t \\ 
     743    F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
     744          &\overline{(\textit{sfx} -        \textit{emp} \lt. S \rt|_{k = 1})}^t 
     745  \end{alignedat} 
    768746\end{equation}  
    769747Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 
     
    783761 
    784762Options are defined through the \ngn{namtra\_qsr} namelist variables. 
    785 When the penetrative solar radiation option is used (\np{ln\_flxqsr}\forcode{ = .true.}), 
     763When the penetrative solar radiation option is used (\np{ln\_flxqsr}~\forcode{= .true.}), 
    786764the solar radiation penetrates the top few tens of meters of the ocean. 
    787 If it is not used (\np{ln\_flxqsr}\forcode{ = .false.}) all the heat flux is absorbed in the first ocean level. 
     765If it is not used (\np{ln\_flxqsr}~\forcode{= .false.}) all the heat flux is absorbed in the first ocean level. 
    788766Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:PE_tra_T} and 
    789767the surface boundary condition is modified to take into account only the non-penetrative part of the surface  
     
    791769\begin{equation} 
    792770  \label{eq:PE_qsr} 
    793   \begin{split} 
    794     \frac{\partial T}{\partial t} &= {\ldots} + \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k}  \\ 
    795     Q_{ns} &= Q_\text{Total} - Q_{sr} 
    796   \end{split} 
     771  \begin{gathered} 
     772    \pd[T]{t} = \ldots + \frac{1}{\rho_o \, C_p \, e_3} \; \pd[I]{k} \\ 
     773    Q_{ns} = Q_\text{Total} - Q_{sr} 
     774  \end{gathered} 
    797775\end{equation} 
    798776where $Q_{sr}$ is the penetrative part of the surface heat flux (\ie the shortwave radiation) and 
    799 $I$ is the downward irradiance ($\left. I \right|_{z=\eta}=Q_{sr}$). 
     777$I$ is the downward irradiance ($\lt. I \rt|_{z = \eta} = Q_{sr}$). 
    800778The additional term in \autoref{eq:PE_qsr} is discretized as follows: 
    801779\begin{equation} 
    802780  \label{eq:tra_qsr} 
    803   \frac{1}{\rho_o\, C_p \,e_3} \; \frac{\partial I}{\partial k} \equiv \frac{1}{\rho_o\, C_p\, e_{3t}} \delta_k \left[ I_w \right] 
     781  \frac{1}{\rho_o \, C_p \, e_3} \, \pd[I]{k} \equiv \frac{1}{\rho_o \, C_p \, e_{3t}} \delta_k [I_w] 
    804782\end{equation} 
    805783 
     
    810788(specified through namelist parameter \np{rn\_abs}). 
    811789It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
    812 of a few tens of centimetres (typically $\xi_0=0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist). 
     790of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn\_si0} in the \ngn{namtra\_qsr} namelist). 
    813791For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to 
    814792larger depths where it contributes to local heating. 
    815793The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen. 
    816 In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}\forcode{ = .true.}) 
     794In the simple 2-waveband light penetration scheme (\np{ln\_qsr\_2bd}~\forcode{= .true.}) 
    817795a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 
    818796leading to the following expression \citep{Paulson1977}: 
    819797\[ 
    820798  % \label{eq:traqsr_iradiance} 
    821   I(z) = Q_{sr} \left[Re^{-z / \xi_0} + \left( 1-R\right) e^{-z / \xi_1} \right] 
     799  I(z) = Q_{sr} \lt[ Re^{- z / \xi_0} + (1 - R) e^{- z / \xi_1} \rt] 
    822800\] 
    823801where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 
    824802It is usually chosen to be 23~m by setting the \np{rn\_si0} namelist parameter. 
    825 The set of default values ($\xi_0$, $\xi_1$, $R$) corresponds to a Type I water in Jerlov's (1968) classification 
     803The set of default values ($\xi_0, \xi_1, R$) corresponds to a Type I water in Jerlov's (1968) classification 
    826804(oligotrophic waters). 
    827805 
     
    840818reproduces quite closely the light penetration profiles predicted by the full spectal model, 
    841819but with much greater computational efficiency. 
    842 The 2-bands formulation does not reproduce the full model very well.  
    843  
    844 The RGB formulation is used when \np{ln\_qsr\_rgb}\forcode{ = .true.}. 
     820The 2-bands formulation does not reproduce the full model very well. 
     821 
     822The RGB formulation is used when \np{ln\_qsr\_rgb}~\forcode{= .true.}. 
    845823The RGB attenuation coefficients (\ie the inverses of the extinction length scales) are tabulated over 
    84682461 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L 
    847825(see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module). 
    848826Four types of chlorophyll can be chosen in the RGB formulation: 
    849 \begin{description}  
    850 \item[\np{nn\_chdta}\forcode{ = 0}] 
     827 
     828\begin{description} 
     829\item[\np{nn\_chdta}~\forcode{= 0}] 
    851830  a constant 0.05 g.Chl/L value everywhere ;  
    852 \item[\np{nn\_chdta}\forcode{ = 1}] 
     831\item[\np{nn\_chdta}~\forcode{= 1}] 
    853832  an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in 
    854833  the vertical direction; 
    855 \item[\np{nn\_chdta}\forcode{ = 2}] 
     834\item[\np{nn\_chdta}~\forcode{= 2}] 
    856835  same as previous case except that a vertical profile of chlorophyl is used. 
    857836  Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value; 
    858 \item[\np{ln\_qsr\_bio}\forcode{ = .true.}] 
     837\item[\np{ln\_qsr\_bio}~\forcode{= .true.}] 
    859838  simulated time varying chlorophyll by TOP biogeochemical model. 
    860839  In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in 
    861   PISCES or LOBSTER and the oceanic heating rate.  
     840  PISCES or LOBSTER and the oceanic heating rate. 
    862841\end{description}  
     842 
    863843The trend in \autoref{eq:tra_qsr} associated with the penetration of the solar radiation is added to 
    864844the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
     
    871851Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor. 
    872852In this case, we have chosen that all remaining radiation is absorbed in the last ocean level 
    873 (\ie $I$ is masked).  
     853(\ie $I$ is masked). 
    874854 
    875855%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    876856\begin{figure}[!t] 
    877857  \begin{center} 
    878     \includegraphics[width=1.0\textwidth]{Fig_TRA_Irradiance} 
     858    \includegraphics[]{Fig_TRA_Irradiance} 
    879859    \caption{ 
    880860      \protect\label{fig:traqsr_irradiance} 
     
    903883\begin{figure}[!t] 
    904884  \begin{center} 
    905     \includegraphics[width=1.0\textwidth]{Fig_TRA_geoth} 
     885    \includegraphics[]{Fig_TRA_geoth} 
    906886    \caption{ 
    907887      \protect\label{fig:geothermal} 
     
    917897This is the default option in \NEMO, and it is implemented using the masking technique. 
    918898However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling. 
    919 This flux is weak compared to surface fluxes (a mean global value of $\sim0.1\;W/m^2$ \citep{Stein_Stein_Nat92}), 
     899This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{Stein_Stein_Nat92}), 
    920900but it warms systematically the ocean and acts on the densest water masses. 
    921901Taking this flux into account in a global ocean model increases the deepest overturning cell 
    922 (\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups  \citep{Emile-Geay_Madec_OS09}.  
     902(\ie the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{Emile-Geay_Madec_OS09}. 
    923903 
    924904Options are defined through the  \ngn{namtra\_bbc} namelist variables. 
     
    939919%-------------------------------------------------------------------------------------------------------------- 
    940920 
    941 Options are defined through the  \ngn{nambbl} namelist variables. 
     921Options are defined through the \ngn{nambbl} namelist variables. 
    942922In a $z$-coordinate configuration, the bottom topography is represented by a series of discrete steps. 
    943923This is not adequate to represent gravity driven downslope flows. 
     
    951931sometimes over a thickness much larger than the thickness of the observed gravity plume. 
    952932A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of 
    953 a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved.  
     933a sill \citep{Willebrand_al_PO01}, and the thickness of the plume is not resolved. 
    954934 
    955935The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{Beckmann_Doscher1997}, 
     
    964944%        Diffusive BBL 
    965945% ------------------------------------------------------------------------------------------------------------- 
    966 \subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}\forcode{ = 1})} 
     946\subsection{Diffusive bottom boundary layer (\protect\np{nn\_bbl\_ldf}~\forcode{= 1})} 
    967947\label{subsec:TRA_bbl_diff} 
    968948 
     
    971951\[ 
    972952  % \label{eq:tra_bbl_diff} 
    973   {\rm {\bf F}}_\sigma=A_l^\sigma \; \nabla_\sigma T 
     953  \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T 
    974954\] 
    975 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 
    976 and  $A_l^\sigma$ the lateral diffusivity in the BBL. 
     955with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and 
     956$A_l^\sigma$ the lateral diffusivity in the BBL. 
    977957Following \citet{Beckmann_Doscher1997}, the latter is prescribed with a spatial dependence, 
    978958\ie in the conditional form 
    979959\begin{equation} 
    980960  \label{eq:tra_bbl_coef} 
    981   A_l^\sigma (i,j,t)=\left\{ { 
    982       \begin{array}{l} 
    983         A_{bbl}  \quad \quad   \mbox{if}  \quad   \nabla_\sigma \rho  \cdot  \nabla H<0 \\ \\ 
    984         0\quad \quad \;\,\mbox{otherwise} \\ 
    985       \end{array}} 
    986   \right. 
    987 \end{equation}  
     961  A_l^\sigma (i,j,t) = 
     962      \begin{cases} 
     963        A_{bbl} & \text{if~} \nabla_\sigma \rho \cdot \nabla H < 0 \\ 
     964        \\ 
     965        0      & \text{otherwise} \\ 
     966      \end{cases} 
     967\end{equation} 
    988968where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn\_ahtbbl} and 
    989969usually set to a value much larger than the one used for lateral mixing in the open ocean. 
     
    995975\[ 
    996976  % \label{eq:tra_bbl_Drho} 
    997   \nabla_\sigma \rho / \rho = \alpha \,\nabla_\sigma T + \beta   \,\nabla_\sigma S 
     977  \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S 
    998978\] 
    999 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline{T}^\sigma$, 
    1000 $\overline{S}^\sigma$ and $\overline{H}^\sigma$, the along bottom mean temperature, salinity and depth, respectively. 
     979where $\rho$, $\alpha$ and $\beta$ are functions of $\overline T^\sigma$, $\overline S^\sigma$ and 
     980$\overline H^\sigma$, the along bottom mean temperature, salinity and depth, respectively. 
    1001981 
    1002982% ------------------------------------------------------------------------------------------------------------- 
    1003983%        Advective BBL 
    1004984% ------------------------------------------------------------------------------------------------------------- 
    1005 \subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}\forcode{ = 1..2})} 
     985\subsection{Advective bottom boundary layer  (\protect\np{nn\_bbl\_adv}~\forcode{= 1..2})} 
    1006986\label{subsec:TRA_bbl_adv} 
    1007987 
    1008 %\sgacomment{"downsloping flow" has been replaced by "downslope flow" in the following 
    1009 %if this is not what is meant then "downwards sloping flow" is also a possibility"} 
     988%\sgacomment{ 
     989%  "downsloping flow" has been replaced by "downslope flow" in the following 
     990%  if this is not what is meant then "downwards sloping flow" is also a possibility" 
     991%} 
    1010992 
    1011993%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1012994\begin{figure}[!t] 
    1013995  \begin{center} 
    1014     \includegraphics[width=0.7\textwidth]{Fig_BBL_adv} 
     996    \includegraphics[]{Fig_BBL_adv} 
    1015997    \caption{ 
    1016998      \protect\label{fig:bbl} 
    1017999      Advective/diffusive Bottom Boundary Layer. 
    1018       The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i+1}_{kdnw}$. 
     1000      The BBL parameterisation is activated when $\rho^i_{kup}$ is larger than $\rho^{i + 1}_{kdnw}$. 
    10191001      Red arrows indicate the additional overturning circulation due to the advective BBL. 
    10201002      The transport of the downslope flow is defined either as the transport of the bottom ocean cell (black arrow), 
     
    10261008%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    10271009 
    1028  
    10291010%!!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity 
    10301011%!!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation 
    1031 %!!        \ie transport proportional to the along-slope density gradient 
     1012%!!        i.e. transport proportional to the along-slope density gradient 
    10321013 
    10331014%%%gmcomment   :  this section has to be really written 
    10341015 
    1035 When applying an advective BBL (\np{nn\_bbl\_adv}\forcode{ = 1..2}), an overturning circulation is added which 
     1016When applying an advective BBL (\np{nn\_bbl\_adv}~\forcode{= 1..2}), an overturning circulation is added which 
    10361017connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope. 
    1037 The density difference causes dense water to move down the slope.  
    1038  
    1039 \np{nn\_bbl\_adv}\forcode{ = 1}: 
     1018The density difference causes dense water to move down the slope. 
     1019 
     1020\np{nn\_bbl\_adv}~\forcode{= 1}: 
    10401021the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step 
    10411022(see black arrow in \autoref{fig:bbl}) \citep{Beckmann_Doscher1997}. 
    10421023It is a \textit{conditional advection}, that is, advection is allowed only 
    1043 if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho  \cdot  \nabla H<0$) and 
    1044 if the velocity is directed towards greater depth (\ie $\vect{U}  \cdot  \nabla H>0$). 
    1045  
    1046 \np{nn\_bbl\_adv}\forcode{ = 2}: 
     1024if dense water overlies less dense water on the slope (\ie $\nabla_\sigma \rho \cdot \nabla H < 0$) and 
     1025if the velocity is directed towards greater depth (\ie $\vect U \cdot \nabla H > 0$). 
     1026 
     1027\np{nn\_bbl\_adv}~\forcode{= 2}: 
    10471028the downslope velocity is chosen to be proportional to $\Delta \rho$, 
    10481029the density difference between the higher cell and lower cell densities \citep{Campin_Goosse_Tel99}. 
    10491030The advection is allowed only  if dense water overlies less dense water on the slope 
    1050 (\ie $\nabla_\sigma \rho  \cdot  \nabla H<0$). 
     1031(\ie $\nabla_\sigma \rho \cdot \nabla H < 0$). 
    10511032For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:bbl}), 
    10521033is simply given by the following expression: 
    10531034\[ 
    10541035  % \label{eq:bbl_Utr} 
    1055   u^{tr}_{bbl} = \gamma \, g \frac{\Delta \rho}{\rho_o}  e_{1u} \; min \left( {e_{3u}}_{kup},{e_{3u}}_{kdwn} \right) 
     1036  u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn}) 
    10561037\] 
    10571038where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as \np{rn\_gambbl}, 
     
    10621043The possible values for $\gamma$ range between 1 and $10~s$ \citep{Campin_Goosse_Tel99}. 
    10631044 
    1064 Scalar properties are advected by this additional transport $( u^{tr}_{bbl}, v^{tr}_{bbl} )$ using the upwind scheme. 
     1045Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme. 
    10651046Such a diffusive advective scheme has been chosen to mimic the entrainment between the downslope plume and 
    10661047the surrounding water at intermediate depths. 
     
    10711052the downslope flow \autoref{eq:bbl_dw}, the horizontal \autoref{eq:bbl_hor} and 
    10721053the upward \autoref{eq:bbl_up} return flows as follows:  
    1073 \begin{align} 
     1054\begin{alignat}{3} 
     1055  \label{eq:bbl_dw} 
    10741056  \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} 
    1075                             +  \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}}  \left( T^{sh}_{kup} - T^{do}_{kdw} \right)  \label{eq:bbl_dw} \\ 
    1076                             % 
     1057                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\ 
     1058  \label{eq:bbl_hor} 
    10771059  \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
    1078                             + \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}}   \left( T^{do}_{kup} - T^{sh}_{kup} \right)   \label{eq:bbl_hor} \\ 
    1079                             % 
     1060                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\ 
     1061  % 
    10801062  \intertext{and for $k =kdw-1,\;..., \; kup$ :} 
    10811063  % 
     1064  \label{eq:bbl_up} 
    10821065  \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k} 
    1083                           + \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   \left( T^{do}_{k+1} - T^{sh}_{k} \right)   \label{eq:bbl_up} 
    1084 \end{align} 
    1085 where $b_t$ is the $T$-cell volume.  
    1086  
    1087 Note that the BBL transport, $( u^{tr}_{bbl}, v^{tr}_{bbl} )$, is available in the model outputs. 
     1066                                &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   &&\lt( T^{do}_{k +1} - T^{sh}_{k}   \rt) 
     1067\end{alignat} 
     1068where $b_t$ is the $T$-cell volume. 
     1069 
     1070Note that the BBL transport, $(u^{tr}_{bbl},v^{tr}_{bbl})$, is available in the model outputs. 
    10881071It has to be used to compute the effective velocity as well as the effective overturning circulation. 
    10891072 
     
    11011084\begin{equation} 
    11021085  \label{eq:tra_dmp} 
    1103   \begin{split} 
    1104     \frac{\partial T}{\partial t}=\;\cdots \;-\gamma \,\left( {T-T_o } \right) \\ 
    1105     \frac{\partial S}{\partial t}=\;\cdots \;-\gamma \,\left( {S-S_o } \right) 
    1106   \end{split} 
     1086  \begin{gathered} 
     1087    \pd[T]{t} = \cdots - \gamma (T - T_o) \\ 
     1088    \pd[S]{t} = \cdots - \gamma (S - S_o) 
     1089  \end{gathered} 
    11071090\end{equation}  
    11081091where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields 
     
    11111094The restoring term is added when the namelist parameter \np{ln\_tradmp} is set to true. 
    11121095It also requires that both \np{ln\_tsd\_init} and \np{ln\_tsd\_tradmp} are set to true in 
    1113 \textit{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set 
     1096\ngn{namtsd} namelist as well as \np{sn\_tem} and \np{sn\_sal} structures are correctly set 
    11141097(\ie that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread}, 
    11151098see \autoref{subsec:SBC_fldread}). 
     
    11281111The second case corresponds to the use of the robust diagnostic method \citep{Sarmiento1982}. 
    11291112It allows us to find the velocity field consistent with the model dynamics whilst 
    1130 having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$).  
     1113having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 
    11311114 
    11321115The robust diagnostic method is very efficient in preventing temperature drift in intermediate waters but 
     
    11401123\citep{Madec_al_JPO96}. 
    11411124 
    1142 \subsection{Generating \ifile{resto} using DMP\_TOOLS} 
    1143  
    1144 DMP\_TOOLS can be used to generate a netcdf file containing the restoration coefficient $\gamma$. 
    1145 Note that in order to maintain bit comparison with previous NEMO versions DMP\_TOOLS must be compiled and 
    1146 run on the same machine as the NEMO model. 
    1147 A \ifile{mesh\_mask} file for the model configuration is required as an input. 
    1148 This can be generated by carrying out a short model run with the namelist parameter \np{nn\_msh} set to 1. 
    1149 The namelist parameter \np{ln\_tradmp} will also need to be set to .false. for this to work. 
    1150 The \ngn{nam\_dmp\_create} namelist in the DMP\_TOOLS directory is used to specify options for 
    1151 the restoration coefficient. 
    1152  
    1153 %--------------------------------------------nam_dmp_create------------------------------------------------- 
    1154 %\namtools{namelist_dmp} 
    1155 %------------------------------------------------------------------------------------------------------- 
    1156  
    1157 \np{cp\_cfg}, \np{cp\_cpz}, \np{jp\_cfg} and \np{jperio} specify the model configuration being used and 
    1158 should be the same as specified in \ngn{namcfg}. 
    1159 The variable \np{lzoom} is used to specify that the damping is being used as in case \textit{a} above to 
    1160 provide boundary conditions to a zoom configuration. 
    1161 In the case of the arctic or antarctic zoom configurations this includes some specific treatment. 
    1162 Otherwise damping is applied to the 6 grid points along the ocean boundaries. 
    1163 The open boundaries are specified by the variables \np{lzoom\_n}, \np{lzoom\_e}, \np{lzoom\_s}, \np{lzoom\_w} in 
    1164 the \ngn{nam\_zoom\_dmp} name list. 
    1165  
    1166 The remaining switch namelist variables determine the spatial variation of the restoration coefficient in 
    1167 non-zoom configurations. 
    1168 \np{ln\_full\_field} specifies that newtonian damping should be applied to the whole model domain. 
    1169 \np{ln\_med\_red\_seas} specifies grid specific restoration coefficients in the Mediterranean Sea for 
    1170 the ORCA4, ORCA2 and ORCA05 configurations. 
    1171 If \np{ln\_old\_31\_lev\_code} is set then the depth variation of the coeffients will be specified as 
    1172 a function of the model number. 
    1173 This option is included to allow backwards compatability of the ORCA2 reference configurations with 
    1174 previous model versions. 
    1175 \np{ln\_coast} specifies that the restoration coefficient should be reduced near to coastlines. 
    1176 This option only has an effect if \np{ln\_full\_field} is true. 
    1177 \np{ln\_zero\_top\_layer} specifies that the restoration coefficient should be zero in the surface layer. 
    1178 Finally \np{ln\_custom} specifies that the custom module will be called. 
    1179 This module is contained in the file \mdl{custom} and can be edited by users. 
    1180 For example damping could be applied in a specific region. 
    1181  
    1182 The restoration coefficient can be set to zero in equatorial regions by 
    1183 specifying a positive value of \np{nn\_hdmp}.  
    1184 Equatorward of this latitude the restoration coefficient will be zero with a smooth transition to  
    1185 the full values of a 10\deg latitud band.  
    1186 This is often used because of the short adjustment time scale in the equatorial region 
    1187 \citep{Reverdin1991, Fujio1991, Marti_PhD92}. 
    1188 The time scale associated with the damping depends on the depth as a hyperbolic tangent, 
    1189 with \np{rn\_surf} as surface value, \np{rn\_bot} as bottom value and a transition depth of \np{rn\_dep}.   
     1125For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under 
     1126\path{./tools/DMP_TOOLS}. 
    11901127 
    11911128% ================================================================ 
     
    11991136%-------------------------------------------------------------------------------------------------------------- 
    12001137 
    1201 Options are defined through the  \ngn{namdom} namelist variables. 
     1138Options are defined through the \ngn{namdom} namelist variables. 
    12021139The general framework for tracer time stepping is a modified leap-frog scheme \citep{Leclair_Madec_OM09}, 
    12031140\ie a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:STP_mLF}): 
    12041141\begin{equation} 
    12051142  \label{eq:tra_nxt} 
    1206   \begin{aligned} 
    1207     (e_{3t}T)^{t+\rdt} &= (e_{3t}T)_f^{t-\rdt} &+ 2 \, \rdt  \,e_{3t}^t\ \text{RHS}^t &   \\ \\ 
    1208     (e_{3t}T)_f^t  \;\ \quad &= (e_{3t}T)^t \;\quad 
    1209     &+\gamma \,\left[ {(e_{3t}T)_f^{t-\rdt} -2(e_{3t}T)^t+(e_{3t}T)^{t+\rdt}} \right] &  \\ 
    1210     & &- \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]  & 
    1211   \end{aligned} 
     1143  \begin{alignedat}{3} 
     1144    &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\ 
     1145    &(e_{3t}T)_f^t        &&= (e_{3t}T)^t            &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\ 
     1146    &                     &&                         &&- \, \gamma \, \rdt \, \lt[ Q^{t + \rdt/2} - Q^{t - \rdt/2} \rt]   
     1147  \end{alignedat} 
    12121148\end{equation}  
    12131149where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values, 
     
    12151151(\ie fluxes plus content in mass exchanges). 
    12161152$\gamma$ is initialized as \np{rn\_atfp} (\textbf{namelist} parameter). 
    1217 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3}. 
     1153Its default value is \np{rn\_atfp}~\forcode{= 10.e-3}. 
    12181154Note that the forcing correction term in the filter is not applied in linear free surface 
    1219 (\jp{lk\_vvl}\forcode{ = .false.}) (see \autoref{subsec:TRA_sbc}. 
     1155(\jp{lk\_vvl}~\forcode{= .false.}) (see \autoref{subsec:TRA_sbc}). 
    12201156Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$. 
    12211157 
    1222 When the vertical mixing is solved implicitly, 
    1223 the update of the \textit{next} tracer fields is done in module \mdl{trazdf}. 
     1158When the vertical mixing is solved implicitly, the update of the \textit{next} tracer fields is done in 
     1159\mdl{trazdf} module. 
    12241160In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} module. 
    12251161 
    12261162In order to prepare for the computation of the \textit{next} time step, a swap of tracer arrays is performed: 
    1227 $T^{t-\rdt} = T^t$ and $T^t = T_f$.  
     1163$T^{t - \rdt} = T^t$ and $T^t = T_f$. 
    12281164 
    12291165% ================================================================ 
     
    12401176%        Equation of State 
    12411177% ------------------------------------------------------------------------------------------------------------- 
    1242 \subsection{Equation of seawater (\protect\np{nn\_eos}\forcode{ = -1..1})} 
     1178\subsection{Equation of seawater (\protect\np{nn\_eos}~\forcode{= -1..1})} 
    12431179\label{subsec:TRA_eos} 
    12441180 
     
    12641200To that purposed, a simplified EOS (S-EOS) inspired by \citet{Vallis06} is also available. 
    12651201 
    1266 In the computer code, a density anomaly, $d_a= \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. 
     1202In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. 
    12671203Called \textit{rau0} in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 
    12681204This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as, 
     
    12701206density in the World Ocean varies by no more than 2$\%$ from that value \citep{Gill1982}. 
    12711207 
    1272 Options are defined through the  \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which 
     1208Options are defined through the \ngn{nameos} namelist variables, and in particular \np{nn\_eos} which 
    12731209controls the EOS used (\forcode{= -1} for TEOS10 ; \forcode{= 0} for EOS-80 ; \forcode{= 1} for S-EOS). 
     1210 
    12741211\begin{description} 
    1275 \item[\np{nn\_eos}\forcode{ = -1}] 
     1212\item[\np{nn\_eos}~\forcode{= -1}] 
    12761213  the polyTEOS10-bsq equation of seawater \citep{Roquet_OM2015} is used. 
    12771214  The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 
     
    12821219  the TEOS-10 rational function approximation for hydrographic data analysis \citep{TEOS10}. 
    12831220  A key point is that conservative state variables are used: 
    1284   Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{C}, notation: $\Theta$). 
     1221  Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \degC, notation: $\Theta$). 
    12851222  The pressure in decibars is approximated by the depth in meters. 
    12861223  With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. 
    1287   It is set to $C_p=3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 
    1288  
     1224  It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{TEOS10}. 
    12891225  Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$. 
    12901226  In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and 
     
    12931229  either computing the air-sea and ice-sea fluxes (forced mode) or 
    12941230  sending the SST field to the atmosphere (coupled mode). 
    1295  
    1296 \item[\np{nn\_eos}\forcode{ = 0}] 
     1231\item[\np{nn\_eos}~\forcode{= 0}] 
    12971232  the polyEOS80-bsq equation of seawater is used. 
    12981233  It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to 
     
    13051240  pressure \citep{UNESCO1983}. 
    13061241  Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which 
    1307   is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value.  
    1308   
    1309 \item[\np{nn\_eos}\forcode{ = 1}] 
     1242  is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 
     1243\item[\np{nn\_eos}~\forcode{= 1}] 
    13101244  a simplified EOS (S-EOS) inspired by \citet{Vallis06} is chosen, 
    13111245  the coefficients of which has been optimized to fit the behavior of TEOS10 
     
    13171251  as well as between \textit{absolute} and \textit{practical} salinity. 
    13181252  S-EOS takes the following expression: 
    1319   \[ 
     1253  \begin{gather*} 
    13201254    % \label{eq:tra_S-EOS} 
    1321     \begin{split} 
    1322       d_a(T,S,z)  =  ( & - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a  \\ 
    1323       & + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a  \\ 
    1324       & - \nu \; T_a \; S_a \;  ) \; / \; \rho_o                     \\ 
    1325       with \ \  T_a = T-10  \; ;  & \;  S_a = S-35  \; ;\;  \rho_o = 1026~Kg/m^3 
    1326     \end{split} 
    1327   \] 
     1255    \begin{alignedat}{2} 
     1256    &d_a(T,S,z) = \frac{1}{\rho_o} \big[ &- a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * &T_a \big. \\ 
     1257    &                                    &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a       \\   
     1258    &                              \big. &- \nu \;                           T_a                  &S_a \big] \\ 
     1259    \end{alignedat} 
     1260    \\ 
     1261    \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3 
     1262  \end{gather*} 
    13281263  where the computer name of the coefficients as well as their standard value are given in \autoref{tab:SEOS}. 
    1329   In fact, when choosing S-EOS, various approximation of EOS can be specified simply by changing the associated coefficients. 
    1330   Setting to zero the two thermobaric coefficients ($\mu_1$, $\mu_2$) remove thermobaric effect from S-EOS. 
    1331   setting to zero the three cabbeling coefficients ($\lambda_1$, $\lambda_2$, $\nu$) remove cabbeling effect from S-EOS. 
     1264  In fact, when choosing S-EOS, various approximation of EOS can be specified simply by 
     1265  changing the associated coefficients. 
     1266  Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$ remove thermobaric effect from S-EOS. 
     1267  setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$ remove cabbeling effect from 
     1268  S-EOS. 
    13321269  Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 
    13331270\end{description} 
    1334  
    13351271 
    13361272%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    13371273\begin{table}[!tb] 
    13381274  \begin{center} 
    1339     \begin{tabular}{|p{26pt}|p{72pt}|p{56pt}|p{136pt}|} 
     1275    \begin{tabular}{|l|l|l|l|} 
    13401276      \hline 
    1341       coeff.   & computer name   & S-EOS     &  description                      \\ \hline 
    1342       $a_0$       & \np{rn\_a0}     & 1.6550 $10^{-1}$ &  linear thermal expansion coeff.    \\ \hline 
    1343       $b_0$       & \np{rn\_b0}     & 7.6554 $10^{-1}$ &  linear haline  expansion coeff.    \\ \hline 
    1344       $\lambda_1$ & \np{rn\_lambda1}& 5.9520 $10^{-2}$ &  cabbeling coeff. in $T^2$          \\ \hline 
    1345       $\lambda_2$ & \np{rn\_lambda2}& 5.4914 $10^{-4}$ &  cabbeling coeff. in $S^2$       \\ \hline 
    1346       $\nu$       & \np{rn\_nu}     & 2.4341 $10^{-3}$ &  cabbeling coeff. in $T \, S$       \\ \hline 
    1347       $\mu_1$     & \np{rn\_mu1}    & 1.4970 $10^{-4}$ &  thermobaric coeff. in T         \\ \hline 
    1348       $\mu_2$     & \np{rn\_mu2}    & 1.1090 $10^{-5}$ &  thermobaric coeff. in S            \\ \hline 
     1277      coeff.      & computer name   & S-EOS           & description                      \\ 
     1278      \hline 
     1279      $a_0$       & \np{rn\_a0}     & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 
     1280      \hline 
     1281      $b_0$       & \np{rn\_b0}     & $7.6554~10^{-1}$ & linear haline  expansion coeff. \\ 
     1282      \hline 
     1283      $\lambda_1$ & \np{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$       \\ 
     1284      \hline 
     1285      $\lambda_2$ & \np{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$       \\ 
     1286      \hline 
     1287      $\nu$       & \np{rn\_nu}     & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$    \\ 
     1288      \hline 
     1289      $\mu_1$     & \np{rn\_mu1}    & $1.4970~10^{-4}$ & thermobaric coeff. in T         \\ 
     1290      \hline 
     1291      $\mu_2$     & \np{rn\_mu2}    & $1.1090~10^{-5}$ & thermobaric coeff. in S         \\ 
     1292      \hline 
    13491293    \end{tabular} 
    13501294    \caption{ 
     
    13521296      Standard value of S-EOS coefficients. 
    13531297    } 
    1354   \end{center} 
     1298\end{center} 
    13551299\end{table} 
    13561300%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    13571301 
    1358  
    13591302% ------------------------------------------------------------------------------------------------------------- 
    13601303%        Brunt-V\"{a}is\"{a}l\"{a} Frequency 
    13611304% ------------------------------------------------------------------------------------------------------------- 
    1362 \subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}\forcode{ = 0..2})} 
     1305\subsection{Brunt-V\"{a}is\"{a}l\"{a} frequency (\protect\np{nn\_eos}~\forcode{= 0..2})} 
    13631306\label{subsec:TRA_bn2} 
    13641307 
    1365 An accurate computation of the ocean stability (\ie of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of 
     1308An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of 
    13661309paramount importance as determine the ocean stratification and is used in several ocean parameterisations 
    13671310(namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion, 
     
    13721315\[ 
    13731316  % \label{eq:tra_bn2} 
    1374   N^2 =  \frac{g}{e_{3w}} \left(   \beta \;\delta_{k+1/2}[S] - \alpha \;\delta_{k+1/2}[T]   \right) 
     1317  N^2 = \frac{g}{e_{3w}} \lt( \beta \; \delta_{k + 1/2}[S] - \alpha \; \delta_{k + 1/2}[T] \rt) 
    13751318\] 
    1376 where $(T,S) = (\Theta, S_A)$ for TEOS10, $= (\theta, S_p)$ for TEOS-80, or $=(T,S)$ for S-EOS, 
    1377 and, $\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
    1378 The coefficients are a polynomial function of temperature, salinity and depth which 
    1379 expression depends on the chosen EOS. 
     1319where $(T,S) = (\Theta,S_A)$ for TEOS10, $(\theta,S_p)$ for TEOS-80, or $(T,S)$ for S-EOS, and, 
     1320$\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
     1321The coefficients are a polynomial function of temperature, salinity and depth which expression depends on 
     1322the chosen EOS. 
    13801323They are computed through \textit{eos\_rab}, a \fortran function that can be found in \mdl{eosbn2}. 
    13811324 
     
    13901333  \label{eq:tra_eos_fzp} 
    13911334  \begin{split} 
    1392     T_f (S,p) = \left( -0.0575 + 1.710523 \;10^{-3} \, \sqrt{S} -  2.154996 \;10^{-4} \,S  \right) \ S    \\ 
    1393     - 7.53\,10^{-3} \ \ p 
    1394   \end{split} 
     1335    &T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\ 
     1336    &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\  
     1337    &\text{and~} d = -7.53~10^{-3} 
     1338    \end{split} 
    13951339\end{equation} 
    13961340 
    13971341\autoref{eq:tra_eos_fzp} is only used to compute the potential freezing point of sea water 
    1398 (\ie referenced to the surface $p=0$), 
     1342(\ie referenced to the surface $p = 0$), 
    13991343thus the pressure dependent terms in \autoref{eq:tra_eos_fzp} (last term) have been dropped. 
    14001344The freezing point is computed through \textit{eos\_fzp}, 
    1401 a \fortran function that can be found in \mdl{eosbn2}.   
    1402  
     1345a \fortran function that can be found in \mdl{eosbn2}. 
    14031346 
    14041347% ------------------------------------------------------------------------------------------------------------- 
     
    14111354% 
    14121355 
    1413  
    14141356% ================================================================ 
    14151357% Horizontal Derivative in zps-coordinate  
     
    14211363I've changed "derivative" to "difference" and "mean" to "average"} 
    14221364 
    1423 With partial cells (\np{ln\_zps}\forcode{ = .true.}) at bottom and top (\np{ln\_isfcav}\forcode{ = .true.}), 
     1365With partial cells (\np{ln\_zps}~\forcode{= .true.}) at bottom and top (\np{ln\_isfcav}~\forcode{= .true.}), 
    14241366in general, tracers in horizontally adjacent cells live at different depths. 
    14251367Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and 
    14261368the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 
    1427 The partial cell properties at the top (\np{ln\_isfcav}\forcode{ = .true.}) are computed in the same way as 
     1369The partial cell properties at the top (\np{ln\_isfcav}~\forcode{= .true.}) are computed in the same way as 
    14281370for the bottom. 
    14291371So, only the bottom interpolation is explained below. 
     
    14321374a linear interpolation in the vertical is used to approximate the deeper tracer as if 
    14331375it actually lived at the depth of the shallower tracer point (\autoref{fig:Partial_step_scheme}). 
    1434 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde{T}$, is: 
     1376For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde T$, is: 
    14351377 
    14361378%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    14371379\begin{figure}[!p] 
    14381380  \begin{center} 
    1439     \includegraphics[width=0.9\textwidth]{Fig_partial_step_scheme} 
     1381    \includegraphics[]{Fig_partial_step_scheme} 
    14401382    \caption{ 
    14411383      \protect\label{fig:Partial_step_scheme} 
    14421384      Discretisation of the horizontal difference and average of tracers in the $z$-partial step coordinate 
    1443       (\protect\np{ln\_zps}\forcode{ = .true.}) in the case $( e3w_k^{i+1} - e3w_k^i  )>0$. 
    1444       A linear interpolation is used to estimate $\widetilde{T}_k^{i+1}$, 
     1385      (\protect\np{ln\_zps}~\forcode{= .true.}) in the case $(e3w_k^{i + 1} - e3w_k^i) > 0$. 
     1386      A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$, 
    14451387      the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 
    1446       The horizontal difference is then given by: $\delta_{i+1/2} T_k=  \widetilde{T}_k^{\,i+1} -T_k^{\,i}$ and 
    1447       the average by: $\overline{T}_k^{\,i+1/2}= ( \widetilde{T}_k^{\,i+1/2} - T_k^{\,i} ) / 2$. 
     1388      The horizontal difference is then given by: $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and 
     1389      the average by: $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$. 
    14481390    } 
    14491391  \end{center} 
     
    14511393%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    14521394\[ 
    1453   \widetilde{T}= \left\{ 
    1454     \begin{aligned} 
    1455       &T^{\,i+1}      -\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right)}{ e_{3w}^{i+1} }\;\delta_k T^{i+1} 
    1456       && \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$   }  \\ \\ 
    1457       &T^{\,i} \ \ \ \,+\frac{ \left( e_{3w}^{i+1} -e_{3w}^i \right) }{e_{3w}^i       }\;\delta_k T^{i+1} 
    1458       && \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
    1459     \end{aligned} 
    1460   \right. 
     1395  \widetilde T = \lt\{ 
     1396    \begin{alignedat}{2} 
     1397      &T^{\, i + 1} &-\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt) }{ e_{3w}^{i + 1} } \; \delta_k T^{i + 1} 
     1398      & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\ \\ 
     1399      &T^{\, i}     &+\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt )}{e_{3w}^i       } \; \delta_k T^{i + 1} 
     1400      & \quad \text{if $e_{3w}^{i + 1} <    e_{3w}^i$} 
     1401    \end{alignedat} 
     1402  \rt. 
    14611403\] 
    14621404and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are:  
    14631405\begin{equation} 
    14641406  \label{eq:zps_hde} 
    1465   \begin{aligned} 
    1466     \delta_{i+1/2} T= 
     1407  \begin{split} 
     1408    \delta_{i + 1/2} T       &= 
    14671409    \begin{cases} 
    1468       \ \ \ \widetilde {T}\quad\ -T^i     & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\ 
    1469       \ \ \ T^{\,i+1}-\widetilde{T}    & \ \ \quad\quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
     1410                                \widetilde T - T^i          & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
     1411                                \\ 
     1412                                T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} <    e_{3w}^i 
    14701413    \end{cases} 
    1471     \\ \\ 
    1472     \overline {T}^{\,i+1/2} \ = 
     1414    \\ 
     1415    \overline T^{\, i + 1/2} &= 
    14731416    \begin{cases} 
    1474       ( \widetilde {T}\ \ \;\,-T^{\,i})    / 2  & \;\ \ \quad\text{if  $\ e_{3w}^{i+1} \geq e_{3w}^i$ } \\ \\ 
    1475       ( T^{\,i+1}-\widetilde{T} ) / 2     & \;\ \ \quad\text{if  $\ e_{3w}^{i+1}    <   e_{3w}^i$   } 
     1417                                (\widetilde T - T^{\, i}   ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
     1418                                \\ 
     1419                                (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} <   e_{3w}^i 
    14761420    \end{cases} 
    1477   \end{aligned} 
     1421  \end{split} 
    14781422\end{equation} 
    14791423 
    14801424The computation of horizontal derivative of tracers as well as of density is performed once for all at 
    14811425each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. 
    1482 It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde{\rho}$, 
     1426It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde \rho$, 
    14831427is not the same as that used for $T$ and $S$. 
    1484 Instead of forming a linear approximation of density, we compute $\widetilde{\rho }$ from the interpolated values of 
     1428Instead of forming a linear approximation of density, we compute $\widetilde \rho$ from the interpolated values of 
    14851429$T$ and $S$, and the pressure at a $u$-point 
    1486 (in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos} ):  
     1430(in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}):  
    14871431\[ 
    14881432  % \label{eq:zps_hde_rho} 
    1489   \widetilde{\rho } = \rho ( {\widetilde{T},\widetilde {S},z_u }) 
    1490   \quad \text{where }\  z_u = \min \left( {z_T^{i+1} ,z_T^i } \right) 
     1433  \widetilde \rho = \rho (\widetilde T,\widetilde S,z_u) \quad \text{where~} z_u = \min \lt( z_T^{i + 1},z_T^i \rt) 
    14911434\] 
    14921435 
Note: See TracChangeset for help on using the changeset viewer.