New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11799 for NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles/chap_model_basics.tex – NEMO

Ignore:
Timestamp:
2019-10-25T16:27:34+02:00 (5 years ago)
Author:
mocavero
Message:

Update the branch to v4.0.1 of the trunk

Location:
NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc

    • Property svn:externals set to
      ^/utils/badges badges
      ^/utils/logos logos
  • NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex

    • Property svn:ignore deleted
  • NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO

    • Property svn:externals set to
      ^/utils/figures/NEMO figures
  • NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles

    • Property svn:ignore
      •  

        old new  
        1 *.aux 
        2 *.bbl 
        3 *.blg 
        4 *.dvi 
        5 *.fdb* 
        6 *.fls 
        7 *.idx 
         1*.ind 
        82*.ilg 
        9 *.ind 
        10 *.log 
        11 *.maf 
        12 *.mtc* 
        13 *.out 
        14 *.pdf 
        15 *.toc 
        16 _minted-* 
  • NEMO/branches/2019/dev_r11470_HPC_12_mpi3/doc/latex/NEMO/subfiles/chap_model_basics.tex

    r11435 r11799  
    1  
    21\documentclass[../main/NEMO_manual]{subfiles} 
    32 
    43\begin{document} 
    54 
    6 % ================================================================ 
    7 % Chapter 1  Model Basics 
    8 % ================================================================ 
    95\chapter{Model Basics} 
    10 \label{chap:PE} 
     6\label{chap:MB} 
     7 
     8\thispagestyle{plain} 
    119 
    1210\chaptertoc 
    1311 
    14 \newpage 
    15  
    16 % ================================================================ 
    17 % Primitive Equations 
    18 % ================================================================ 
     12\paragraph{Changes record} ~\\ 
     13 
     14{\footnotesize 
     15  \begin{tabular}{l||l|l} 
     16    Release          & Author(s)                                   & Modifications       \\ 
     17    \hline 
     18    {\em        4.0} & {\em Mike Bell                            } & {\em Review       } \\ 
     19    {\em        3.6} & {\em Tim Graham and Gurvan Madec          } & {\em Updates      } \\ 
     20    {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ 
     21  \end{tabular} 
     22} 
     23 
     24\clearpage 
     25 
     26%% ================================================================================================= 
    1927\section{Primitive equations} 
    20 \label{sec:PE_PE} 
    21  
    22 % ------------------------------------------------------------------------------------------------------------- 
    23 %        Vector Invariant Formulation 
    24 % ------------------------------------------------------------------------------------------------------------- 
    25  
     28\label{sec:MB_PE} 
     29 
     30%% ================================================================================================= 
    2631\subsection{Vector invariant formulation} 
    27 \label{subsec:PE_Vector} 
     32\label{subsec:MB_PE_vector} 
    2833 
    2934The ocean is a fluid that can be described to a good approximation by the primitive equations, 
     
    3237plus the following additional assumptions made from scale considerations: 
    3338 
    34 \begin{enumerate} 
    35 \item 
    36   \textit{spherical Earth approximation}: the geopotential surfaces are assumed to be oblate spheriods 
    37   that follow the Earth's bulge; these spheroids are approximated by spheres with 
    38   gravity locally vertical (parallel to the Earth's radius) and independent of latitude 
     39\begin{labeling}{Neglect of additional Coriolis terms} 
     40\item [\textit{Spherical Earth approximation}] The geopotential surfaces are assumed to 
     41  be oblate spheroids that follow the Earth's bulge; 
     42  these spheroids are approximated by spheres with gravity locally vertical 
     43  (parallel to the Earth's radius) and independent of latitude 
    3944  \citep[][section 2]{white.hoskins.ea_QJRMS05}. 
    40 \item 
    41   \textit{thin-shell approximation}: the ocean depth is neglected compared to the earth's radius 
    42 \item 
    43   \textit{turbulent closure hypothesis}: the turbulent fluxes 
     45\item [\textit{Thin-shell approximation}] The ocean depth is neglected compared to the earth's radius 
     46\item [\textit{Turbulent closure hypothesis}] The turbulent fluxes 
    4447  (which represent the effect of small scale processes on the large-scale) 
    4548  are expressed in terms of large-scale features 
    46 \item 
    47   \textit{Boussinesq hypothesis}: density variations are neglected except in their contribution to 
    48   the buoyancy force 
     49\item [\textit{Boussinesq hypothesis}] Density variations are neglected except in 
     50  their contribution to the buoyancy force 
    4951  \begin{equation} 
    50     \label{eq:PE_eos} 
     52    \label{eq:MB_PE_eos} 
    5153    \rho = \rho \ (T,S,p) 
    5254  \end{equation} 
    53 \item 
    54   \textit{Hydrostatic hypothesis}: the vertical momentum equation is reduced to a balance between 
    55   the vertical pressure gradient and the buoyancy force 
     55\item [\textit{Hydrostatic hypothesis}] The vertical momentum equation is reduced to 
     56  a balance between the vertical pressure gradient and the buoyancy force 
    5657  (this removes convective processes from the initial Navier-Stokes equations and so 
    5758  convective processes must be parameterized instead) 
    5859  \begin{equation} 
    59     \label{eq:PE_hydrostatic} 
     60    \label{eq:MB_PE_hydrostatic} 
    6061    \pd[p]{z} = - \rho \ g 
    6162  \end{equation} 
    62 \item 
    63   \textit{Incompressibility hypothesis}: the three dimensional divergence of the velocity vector $\vect U$ 
    64   is assumed to be zero. 
     63\item [\textit{Incompressibility hypothesis}] The three dimensional divergence of 
     64  the velocity vector $\vect U$ is assumed to be zero. 
    6565  \begin{equation} 
    66     \label{eq:PE_continuity} 
     66    \label{eq:MB_PE_continuity} 
    6767    \nabla \cdot \vect U = 0 
    6868  \end{equation} 
    69  \item 
    70   \textit{Neglect of additional Coriolis terms}: the Coriolis terms that vary with the cosine of latitude are neglected. 
    71   These terms may be non-negligible where the Brunt-Vaisala frequency $N$ is small, either in the deep ocean or 
    72   in the sub-mesoscale motions of the mixed layer, or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. 
    73   They can be consistently included as part of the ocean dynamics \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are 
    74   retained in the MIT ocean model. 
    75 \end{enumerate} 
     69\item [\textit{Neglect of additional Coriolis terms}] The Coriolis terms that vary with 
     70  the cosine of latitude are neglected. 
     71  These terms may be non-negligible where the Brunt-V\"{a}is\"{a}l\"{a} frequency $N$ is small, 
     72  either in the deep ocean or in the sub-mesoscale motions of the mixed layer, 
     73  or near the equator \citep[][section 1]{white.hoskins.ea_QJRMS05}. 
     74  They can be consistently included as part of the ocean dynamics 
     75  \citep[][section 3(d)]{white.hoskins.ea_QJRMS05} and are retained in the MIT ocean model. 
     76\end{labeling} 
    7677 
    7778Because the gravitational force is so dominant in the equations of large-scale motions, 
     
    7980$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 
    8081\ie\ tangent to the geopotential surfaces. 
    81 Let us define the following variables: $\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 
     82Let us define the following variables: 
     83$\vect U$ the vector velocity, $\vect U = \vect U_h + w \, \vect k$ 
    8284(the subscript $h$ denotes the local horizontal vector, \ie\ over the $(i,j)$ plane), 
    8385$T$ the potential temperature, $S$ the salinity, $\rho$ the \textit{in situ} density. 
     
    8587the following equations: 
    8688\begin{subequations} 
    87   \label{eq:PE} 
     89  \label{eq:MB_PE} 
    8890  \begin{gather} 
    89     \intertext{$-$ the momentum balance} 
    90     \label{eq:PE_dyn} 
    91     \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h 
    92                         - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p 
    93                         + \vect D^{\vect U} + \vect F^{\vect U} \\ 
    94     \intertext{$-$ the heat and salt conservation equations} 
    95     \label{eq:PE_tra_T} 
     91    \shortintertext{$-$ the momentum balance} 
     92    \label{eq:MB_PE_dyn} 
     93    \pd[\vect U_h]{t} = - \lt[ (\nabla \times \vect U) \times \vect U + \frac{1}{2} \nabla \lt( \vect U^2 \rt) \rt]_h - f \; k \times \vect U_h - \frac{1}{\rho_o} \nabla_h p + \vect D^{\vect U} + \vect F^{\vect U} 
     94    \shortintertext{$-$ the heat and salt conservation equations} 
     95    \label{eq:MB_PE_tra_T} 
    9696    \pd[T]{t} = - \nabla \cdot (T \ \vect U) + D^T + F^T \\ 
    97     \label{eq:PE_tra_S} 
     97    \label{eq:MB_PE_tra_S} 
    9898    \pd[S]{t} = - \nabla \cdot (S \ \vect U) + D^S + F^S 
    9999  \end{gather} 
     
    101101where $\nabla$ is the generalised derivative vector operator in $(i,j,k)$ directions, $t$ is the time, 
    102102$z$ is the vertical coordinate, $\rho$ is the \textit{in situ} density given by the equation of state 
    103 (\autoref{eq:PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 
     103(\autoref{eq:MB_PE_eos}), $\rho_o$ is a reference density, $p$ the pressure, 
    104104$f = 2 \vect \Omega \cdot k$ is the Coriolis acceleration 
    105 (where $\vect \Omega$ is the Earth's angular velocity vector), and $g$ is the gravitational acceleration. 
     105(where $\vect \Omega$ is the Earth's angular velocity vector), 
     106and $g$ is the gravitational acceleration. 
    106107$\vect D^{\vect U}$, $D^T$ and $D^S$ are the parameterisations of small-scale physics for momentum, 
    107108temperature and salinity, and $\vect F^{\vect U}$, $F^T$ and $F^S$ surface forcing terms. 
    108 Their nature and formulation are discussed in \autoref{sec:PE_zdf_ldf} and \autoref{subsec:PE_boundary_condition}. 
    109  
    110 % ------------------------------------------------------------------------------------------------------------- 
    111 % Boundary condition 
    112 % ------------------------------------------------------------------------------------------------------------- 
     109Their nature and formulation are discussed in \autoref{sec:MB_zdf_ldf} and 
     110\autoref{subsec:MB_boundary_condition}. 
     111 
     112%% ================================================================================================= 
    113113\subsection{Boundary conditions} 
    114 \label{subsec:PE_boundary_condition} 
     114\label{subsec:MB_boundary_condition} 
    115115 
    116116An ocean is bounded by complex coastlines, bottom topography at its base and 
     
    119119where $H$ is the depth of the ocean bottom and $\eta$ is the height of the sea surface 
    120120(discretisation can introduce additional artificial ``side-wall'' boundaries). 
    121 Both $H$ and $\eta$ are referenced to a surface of constant geopotential (\ie\ a mean sea surface height) on which $z = 0$. 
    122 (\autoref{fig:ocean_bc}). 
    123 Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, and momentum with 
    124 the solid earth, the continental margins, the sea ice and the atmosphere. 
    125 However, some of these fluxes are so weak that even on climatic time scales of thousands of years 
    126 they can be neglected. 
     121Both $H$ and $\eta$ are referenced to a surface of constant geopotential 
     122(\ie\ a mean sea surface height) on which $z = 0$ (\autoref{fig:MB_ocean_bc}). 
     123Through these two boundaries, the ocean can exchange fluxes of heat, fresh water, salt, 
     124and momentum with the solid earth, the continental margins, the sea ice and the atmosphere. 
     125However, some of these fluxes are so weak that 
     126even on climatic time scales of thousands of years they can be neglected. 
    127127In the following, we briefly review the fluxes exchanged at the interfaces between the ocean and 
    128128the other components of the earth system. 
    129129 
    130 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    131 \begin{figure}[!ht] 
    132   \begin{center} 
    133     \includegraphics[width=\textwidth]{Fig_I_ocean_bc} 
    134     \caption{ 
    135       \protect\label{fig:ocean_bc} 
    136       The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 
    137       where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 
    138       Both $H$ and $\eta$ are referenced to $z = 0$. 
    139     } 
    140   \end{center} 
     130\begin{figure} 
     131  \centering 
     132  \includegraphics[width=0.66\textwidth]{Fig_I_ocean_bc} 
     133  \caption[Ocean boundary conditions]{ 
     134    The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 
     135    where $H$ is the depth of the sea floor and $\eta$ the height of the sea surface. 
     136    Both $H$ and $\eta$ are referenced to $z = 0$.} 
     137  \label{fig:MB_ocean_bc} 
    141138\end{figure} 
    142 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    143139 
    144140\begin{description} 
    145 \item[Land - ocean interface:] 
    146   the major flux between continental margins and the ocean is a mass exchange of fresh water through river runoff. 
     141\item [Land - ocean] The major flux between continental margins and the ocean is a mass exchange of 
     142  fresh water through river runoff. 
    147143  Such an exchange modifies the sea surface salinity especially in the vicinity of major river mouths. 
    148   It can be neglected for short range integrations but has to be taken into account for long term integrations as 
     144  It can be neglected for short range integrations but 
     145  has to be taken into account for long term integrations as 
    149146  it influences the characteristics of water masses formed (especially at high latitudes). 
    150147  It is required in order to close the water cycle of the climate system. 
    151   It is usually specified as a fresh water flux at the air-sea interface in the vicinity of river mouths. 
    152 \item[Solid earth - ocean interface:] 
    153   heat and salt fluxes through the sea floor are small, except in special areas of little extent. 
    154   They are usually neglected in the model 
    155   \footnote{ 
     148  It is usually specified as a fresh water flux at the air-sea interface in 
     149  the vicinity of river mouths. 
     150\item [Solid earth - ocean] Heat and salt fluxes through the sea floor are small, 
     151  except in special areas of little extent. 
     152  They are usually neglected in the model \footnote{ 
    156153    In fact, it has been shown that the heat flux associated with the solid Earth cooling 
    157     (\ie\ the geothermal heating) is not negligible for the thermohaline circulation of the world ocean 
    158     (see \autoref{subsec:TRA_bbc}). 
     154    (\ie\ the geothermal heating) is not negligible for 
     155    the thermohaline circulation of the world ocean (see \autoref{subsec:TRA_bbc}). 
    159156  }. 
    160157  The boundary condition is thus set to no flux of heat and salt across solid boundaries. 
    161   For momentum, the situation is different. There is no flow across solid boundaries, 
    162   \ie\ the velocity normal to the ocean bottom and coastlines is zero (in other words, 
    163   the bottom velocity is parallel to solid boundaries). This kinematic boundary condition 
    164   can be expressed as: 
     158  For momentum, the situation is different. 
     159  There is no flow across solid boundaries, 
     160  \ie\ the velocity normal to the ocean bottom and coastlines is zero 
     161  (in other words, the bottom velocity is parallel to solid boundaries). 
     162  This kinematic boundary condition can be expressed as: 
    165163  \begin{equation} 
    166     \label{eq:PE_w_bbc} 
     164    \label{eq:MB_w_bbc} 
    167165    w = - \vect U_h \cdot \nabla_h (H) 
    168166  \end{equation} 
    169167  In addition, the ocean exchanges momentum with the earth through frictional processes. 
    170168  Such momentum transfer occurs at small scales in a boundary layer. 
    171   It must be parameterized in terms of turbulent fluxes using bottom and/or lateral boundary conditions. 
     169  It must be parameterized in terms of turbulent fluxes using 
     170  bottom and/or lateral boundary conditions. 
    172171  Its specification depends on the nature of the physical parameterisation used for 
    173   $\vect D^{\vect U}$ in \autoref{eq:PE_dyn}. 
    174   It is discussed in \autoref{eq:PE_zdf}.% and Chap. III.6 to 9. 
    175 \item[Atmosphere - ocean interface:] 
    176   the kinematic surface condition plus the mass flux of fresh water PE (the precipitation minus evaporation budget) 
    177   leads to: 
     172  $\vect D^{\vect U}$ in \autoref{eq:MB_PE_dyn}. 
     173  It is discussed in \autoref{eq:MB_zdf}. % and Chap. III.6 to 9. 
     174\item [Atmosphere - ocean] The kinematic surface condition plus the mass flux of fresh water PE 
     175  (the precipitation minus evaporation budget) leads to: 
    178176  \[ 
    179     % \label{eq:PE_w_sbc} 
     177    % \label{eq:MB_w_sbc} 
    180178    w = \pd[\eta]{t} + \lt. \vect U_h \rt|_{z = \eta} \cdot \nabla_h (\eta) + P - E 
    181179  \] 
    182   The dynamic boundary condition, neglecting the surface tension (which removes capillary waves from the system) 
    183   leads to the continuity of pressure across the interface $z = \eta$. 
     180  The dynamic boundary condition, neglecting the surface tension 
     181  (which removes capillary waves from the system) leads to 
     182  the continuity of pressure across the interface $z = \eta$. 
    184183  The atmosphere and ocean also exchange horizontal momentum (wind stress), and heat. 
    185 \item[Sea ice - ocean interface:] 
    186   the ocean and sea ice exchange heat, salt, fresh water and momentum. 
     184\item [Sea ice - ocean] The ocean and sea ice exchange heat, salt, fresh water and momentum. 
    187185  The sea surface temperature is constrained to be at the freezing point at the interface. 
    188186  Sea ice salinity is very low ($\sim4-6 \, psu$) compared to those of the ocean ($\sim34 \, psu$). 
    189   The cycle of freezing/melting is associated with fresh water and salt fluxes that cannot be neglected. 
     187  The cycle of freezing/melting is associated with fresh water and salt fluxes that 
     188  cannot be neglected. 
    190189\end{description} 
    191190 
    192 % ================================================================ 
    193 % The Horizontal Pressure Gradient 
    194 % ================================================================ 
     191%% ================================================================================================= 
    195192\section{Horizontal pressure gradient} 
    196 \label{sec:PE_hor_pg} 
    197  
    198 % ------------------------------------------------------------------------------------------------------------- 
    199 % Pressure Formulation 
    200 % ------------------------------------------------------------------------------------------------------------- 
     193\label{sec:MB_hor_pg} 
     194 
     195%% ================================================================================================= 
    201196\subsection{Pressure formulation} 
    202 \label{subsec:PE_p_formulation} 
     197\label{subsec:MB_p_formulation} 
    203198 
    204199The total pressure at a given depth $z$ is composed of a surface pressure $p_s$ at 
    205200a reference geopotential surface ($z = 0$) and a hydrostatic pressure $p_h$ such that: 
    206201$p(i,j,k,t) = p_s(i,j,t) + p_h(i,j,k,t)$. 
    207 The latter is computed by integrating (\autoref{eq:PE_hydrostatic}), 
    208 assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:PE_eos}). 
     202The latter is computed by integrating (\autoref{eq:MB_PE_hydrostatic}), 
     203assuming that pressure in decibars can be approximated by depth in meters in (\autoref{eq:MB_PE_eos}). 
    209204The hydrostatic pressure is then given by: 
    210205\[ 
    211   % \label{eq:PE_pressure} 
     206  % \label{eq:MB_pressure} 
    212207  p_h (i,j,z,t) = \int_{\varsigma = z}^{\varsigma = 0} g \; \rho (T,S,\varsigma) \; d \varsigma 
    213208\] 
    214209Two strategies can be considered for the surface pressure term: 
    215 $(a)$ introduce of a  new variable $\eta$, the free-surface elevation, 
     210\begin{enumerate*}[label=(\textit{\alph*})] 
     211\item introduce of a new variable $\eta$, the free-surface elevation, 
    216212for which a prognostic equation can be established and solved; 
    217 $(b)$ assume that the ocean surface is a rigid lid, 
     213\item assume that the ocean surface is a rigid lid, 
    218214on which the pressure (or its horizontal gradient) can be diagnosed. 
     215\end{enumerate*} 
    219216When the former strategy is used, one solution of the free-surface elevation consists of 
    220217the excitation of external gravity waves. 
     
    227224modifies certain other longwave dynamics (\eg\ barotropic Rossby or planetary waves). 
    228225The rigid-lid hypothesis is an obsolescent feature in modern OGCMs. 
    229 It has been available until the release 3.1 of \NEMO, and it has been removed in release 3.2 and followings. 
     226It has been available until the release 3.1 of \NEMO, 
     227and it has been removed in release 3.2 and followings. 
    230228Only the free surface formulation is now described in this document (see the next sub-section). 
    231229 
    232 % ------------------------------------------------------------------------------------------------------------- 
    233 % Free Surface Formulation 
    234 % ------------------------------------------------------------------------------------------------------------- 
     230%% ================================================================================================= 
    235231\subsection{Free surface formulation} 
    236 \label{subsec:PE_free_surface} 
     232\label{subsec:MB_free_surface} 
    237233 
    238234In the free surface formulation, a variable $\eta$, the sea-surface height, 
    239235is introduced which describes the shape of the air-sea interface. 
    240 This variable is solution of a prognostic equation which is established by forming the vertical average of 
    241 the kinematic surface condition (\autoref{eq:PE_w_bbc}): 
     236This variable is solution of a prognostic equation which is established by 
     237forming the vertical average of the kinematic surface condition (\autoref{eq:MB_w_bbc}): 
    242238\begin{equation} 
    243   \label{eq:PE_ssh} 
     239  \label{eq:MB_ssh} 
    244240  \pd[\eta]{t} = - D + P - E \quad \text{where} \quad D = \nabla \cdot \lt[ (H + \eta) \; \overline{U}_h \, \rt] 
    245241\end{equation} 
    246 and using (\autoref{eq:PE_hydrostatic}) the surface pressure is given by: $p_s = \rho \, g \, \eta$. 
    247  
    248 Allowing the air-sea interface to move introduces the external gravity waves (EGWs) as 
     242and using (\autoref{eq:MB_PE_hydrostatic}) the surface pressure is given by: 
     243$p_s = \rho \, g \, \eta$. 
     244 
     245Allowing the air-sea interface to move introduces 
     246the \textbf{E}xternal \textbf{G}ravity \textbf{W}aves (EGWs) as 
    249247a class of solution of the primitive equations. 
    250248These waves are barotropic (\ie\ nearly independent of depth) and their phase speed is quite high. 
     
    253251Two choices can be made regarding the implementation of the free surface in the model, 
    254252depending on the physical processes of interest. 
    255  
    256 $\bullet$ If one is interested in EGWs, in particular the tides and their interaction with 
    257 the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 
    258 then a non linear free surface is the most appropriate. 
    259 This means that no approximation is made in \autoref{eq:PE_ssh} and that 
    260 the variation of the ocean volume is fully taken into account. 
    261 Note that in order to study the fast time scales associated with EGWs it is necessary to 
    262 minimize time filtering effects 
    263 (use an explicit time scheme with very small time step, or a split-explicit scheme with reasonably small time step, 
    264 see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}). 
    265  
    266 $\bullet$ If one is not interested in EGW but rather sees them as high frequency noise, 
    267 it is possible to apply an explicit filter to slow down the fastest waves while 
    268 not altering the slow barotropic Rossby waves. 
    269 If further, an approximative conservation of heat and salt contents is sufficient for the problem solved, 
    270 then it is sufficient to solve a linearized version of \autoref{eq:PE_ssh}, 
    271 which still allows to take into account freshwater fluxes applied at the ocean surface \citep{roullet.madec_JGR00}. 
    272 Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 
    273  
    274 The filtering of EGWs in models with a free surface is usually a matter of discretisation of 
    275 the temporal derivatives, 
     253\begin{itemize} 
     254\item If one is interested in EGWs, in particular the tides and their interaction with 
     255  the baroclinic structure of the ocean (internal waves) possibly in shallow seas, 
     256  then a non linear free surface is the most appropriate. 
     257  This means that no approximation is made in \autoref{eq:MB_ssh} and that 
     258  the variation of the ocean volume is fully taken into account. 
     259  Note that in order to study the fast time scales associated with EGWs it is necessary to 
     260  minimize time filtering effects 
     261  (use an explicit time scheme with very small time step, 
     262  or a split-explicit scheme with reasonably small time step, 
     263  see \autoref{subsec:DYN_spg_exp} or \autoref{subsec:DYN_spg_ts}). 
     264\item If one is not interested in EGWs but rather sees them as high frequency noise, 
     265  it is possible to apply an explicit filter to slow down the fastest waves while 
     266  not altering the slow barotropic Rossby waves. 
     267  If further, an approximative conservation of heat and salt contents is sufficient for 
     268  the problem solved, then it is sufficient to solve a linearized version of \autoref{eq:MB_ssh}, 
     269  which still allows to take into account freshwater fluxes applied at the ocean surface 
     270  \citep{roullet.madec_JGR00}. 
     271  Nevertheless, with the linearization, an exact conservation of heat and salt contents is lost. 
     272\end{itemize} 
     273The filtering of EGWs in models with a free surface is usually 
     274a matter of discretisation of the temporal derivatives, 
    276275using a split-explicit method \citep{killworth.webb.ea_JPO91, zhang.endoh_JGR92} or 
    277 the implicit scheme \citep{dukowicz.smith_JGR94} or 
    278 the addition of a filtering force in the momentum equation \citep{roullet.madec_JGR00}. 
    279 With the present release, \NEMO\  offers the choice between 
    280 an explicit free surface (see \autoref{subsec:DYN_spg_exp}) or 
    281 a split-explicit scheme strongly inspired the one proposed by \citet{shchepetkin.mcwilliams_OM05} 
    282 (see \autoref{subsec:DYN_spg_ts}). 
    283  
    284 % ================================================================ 
    285 % Curvilinear z-coordinate System 
    286 % ================================================================ 
    287 \section{Curvilinear \textit{z-}coordinate system} 
    288 \label{sec:PE_zco} 
    289  
    290 % ------------------------------------------------------------------------------------------------------------- 
    291 % Tensorial Formalism 
    292 % ------------------------------------------------------------------------------------------------------------- 
     276the implicit scheme \citep{dukowicz.smith_JGR94} or the addition of a filtering force in 
     277the momentum equation \citep{roullet.madec_JGR00}. 
     278With the present release, \NEMO\ offers the choice between an explicit free surface 
     279(see \autoref{subsec:DYN_spg_exp}) or a split-explicit scheme strongly inspired the one proposed by 
     280\citet{shchepetkin.mcwilliams_OM05} (see \autoref{subsec:DYN_spg_ts}). 
     281 
     282%% ================================================================================================= 
     283\section{Curvilinear \textit{z}-coordinate system} 
     284\label{sec:MB_zco} 
     285 
     286%% ================================================================================================= 
    293287\subsection{Tensorial formalism} 
    294 \label{subsec:PE_tensorial} 
     288\label{subsec:MB_tensorial} 
    295289 
    296290In many ocean circulation problems, the flow field has regions of enhanced dynamics 
     
    304298A solution consists of introducing an appropriate coordinate transformation that 
    305299shifts the singular point onto land \citep{madec.imbard_CD96, murray_JCP96}. 
    306 As a consequence, it is important to solve the primitive equations in various curvilinear coordinate systems. 
    307 An efficient way of introducing an appropriate coordinate transform can be found when using a tensorial formalism. 
     300As a consequence, 
     301it is important to solve the primitive equations in various curvilinear coordinate systems. 
     302An efficient way of introducing an appropriate coordinate transform can be found when 
     303using a tensorial formalism. 
    308304This formalism is suited to any multidimensional curvilinear coordinate system. 
    309 Ocean modellers mainly use three-dimensional orthogonal grids on the sphere (spherical earth approximation), 
    310 with preservation of the local vertical. Here we give the simplified equations for this particular case. 
    311 The general case is detailed by \citet{eiseman.stone_SR80} in their survey of the conservation laws of fluid dynamics. 
     305Ocean modellers mainly use three-dimensional orthogonal grids on the sphere 
     306(spherical earth approximation), with preservation of the local vertical. 
     307Here we give the simplified equations for this particular case. 
     308The general case is detailed by \citet{eiseman.stone_SR80} in 
     309their survey of the conservation laws of fluid dynamics. 
    312310 
    313311Let $(i,j,k)$ be a set of orthogonal curvilinear coordinates on 
     
    315313$(i,j,k)$ linked to the earth such that 
    316314$k$ is the local upward vector and $(i,j)$ are two vectors orthogonal to $k$, 
    317 \ie\ along geopotential surfaces (\autoref{fig:referential}). 
     315\ie\ along geopotential surfaces (\autoref{fig:MB_referential}). 
    318316Let $(\lambda,\varphi,z)$ be the geographical coordinate system in which a position is defined by 
    319317the latitude $\varphi(i,j)$, the longitude $\lambda(i,j)$ and 
    320318the distance from the centre of the earth $a + z(k)$ where $a$ is the earth's radius and 
    321 $z$ the altitude above a reference sea level (\autoref{fig:referential}). 
     319$z$ the altitude above a reference sea level (\autoref{fig:MB_referential}). 
    322320The local deformation of the curvilinear coordinate system is given by $e_1$, $e_2$ and $e_3$, 
    323321the three scale factors: 
    324322\begin{equation} 
    325   \label{eq:scale_factors} 
    326   \begin{aligned} 
    327     e_1 &= (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \\ 
    328     e_2 &= (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \\ 
    329     e_3 &= \lt( \pd[z]{k} \rt) 
    330   \end{aligned} 
     323  \label{eq:MB_scale_factors} 
     324    e_1 = (a + z) \lt[ \lt( \pd[\lambda]{i} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{i} \rt)^2 \rt]^{1/2} \quad e_2 = (a + z) \lt[ \lt( \pd[\lambda]{j} \cos \varphi \rt)^2 + \lt( \pd[\varphi]{j} \rt)^2 \rt]^{1/2} \quad e_3 = \lt( \pd[z]{k} \rt) 
    331325\end{equation} 
    332326 
    333 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    334 \begin{figure}[!tb] 
    335   \begin{center} 
    336     \includegraphics[width=\textwidth]{Fig_I_earth_referential} 
    337     \caption{ 
    338       \protect\label{fig:referential} 
    339       the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
    340       coordinate system $(i,j,k)$. 
    341     } 
    342   \end{center} 
     327\begin{figure} 
     328  \centering 
     329  \includegraphics[width=0.33\textwidth]{Fig_I_earth_referential} 
     330  \caption[Geographical and curvilinear coordinate systems]{ 
     331    the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
     332    coordinate system $(i,j,k)$.} 
     333  \label{fig:MB_referential} 
    343334\end{figure} 
    344 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    345335 
    346336Since the ocean depth is far smaller than the earth's radius, $a + z$, can be replaced by $a$ in 
    347 (\autoref{eq:scale_factors}) (thin-shell approximation). 
     337(\autoref{eq:MB_scale_factors}) (thin-shell approximation). 
    348338The resulting horizontal scale factors $e_1$, $e_2$  are independent of $k$ while 
    349339the vertical scale factor is a single function of $k$ as $k$ is parallel to $z$. 
    350340The scalar and vector operators that appear in the primitive equations 
    351 (\autoref{eq:PE_dyn} to \autoref{eq:PE_eos}) can then be written in the tensorial form, 
     341(\autoref{eq:MB_PE_dyn} to \autoref{eq:MB_PE_eos}) can then be written in the tensorial form, 
    352342invariant in any orthogonal horizontal curvilinear coordinate system transformation: 
    353343\begin{subequations} 
    354   % \label{eq:PE_discrete_operators} 
    355   \begin{gather} 
    356     \label{eq:PE_grad} 
    357     \nabla q =   \frac{1}{e_1} \pd[q]{i} \; \vect i 
    358                + \frac{1}{e_2} \pd[q]{j} \; \vect j 
    359                + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 
    360     \label{eq:PE_div} 
    361     \nabla \cdot \vect A =   \frac{1}{e_1 \; e_2} \lt[ \pd[(e_2 \; a_1)]{\partial i} + \pd[(e_1 \; a_2)]{j} \rt] 
    362                            + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] 
    363   \end{gather} 
    364   \begin{multline} 
    365     \label{eq:PE_curl} 
    366       \nabla \times \vect{A} =   \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k}   \rt] \vect i \\ 
    367                                + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i}   \rt] \vect j \\ 
    368                                + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k 
    369   \end{multline} 
    370   \begin{gather} 
    371     \label{eq:PE_lap} 
    372     \Delta q = \nabla \cdot (\nabla q) \\ 
    373     \label{eq:PE_lap_vector} 
    374     \Delta \vect A = \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) 
    375   \end{gather} 
     344  % \label{eq:MB_discrete_operators} 
     345  \begin{align} 
     346    \label{eq:MB_grad} 
     347    \nabla q &=   \frac{1}{e_1} \pd[q]{i} \; \vect i + \frac{1}{e_2} \pd[q]{j} \; \vect j + \frac{1}{e_3} \pd[q]{k} \; \vect k \\ 
     348    \label{eq:MB_div} 
     349    \nabla \cdot \vect A &=   \frac{1}{e_1 \; e_2} \lt[ \pd[(e_2 \; a_1)]{\partial i} + \pd[(e_1 \; a_2)]{j} \rt] + \frac{1}{e_3} \lt[ \pd[a_3]{k} \rt] \\ 
     350    \label{eq:MB_curl} 
     351      \nabla \times \vect{A} &=   \lt[ \frac{1}{e_2} \pd[a_3]{j} - \frac{1}{e_3} \pd[a_2]{k}   \rt] \vect i + \lt[ \frac{1}{e_3} \pd[a_1]{k} - \frac{1}{e_1} \pd[a_3]{i}   \rt] \vect j + \frac{1}{e_1 e_2} \lt[ \pd[(e_2 a_2)]{i} - \pd[(e_1 a_1)]{j} \rt] \vect k \\ 
     352    \label{eq:MB_lap} 
     353    \Delta q &= \nabla \cdot (\nabla q) \\ 
     354    \label{eq:MB_lap_vector} 
     355    \Delta \vect A &= \nabla (\nabla \cdot \vect A) - \nabla \times (\nabla \times \vect A) 
     356  \end{align} 
    376357\end{subequations} 
    377 where $q$ is a scalar quantity and $\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. 
    378  
    379 % ------------------------------------------------------------------------------------------------------------- 
    380 % Continuous Model Equations 
    381 % ------------------------------------------------------------------------------------------------------------- 
     358where $q$ is a scalar quantity and 
     359$\vect A = (a_1,a_2,a_3)$ a vector in the $(i,j,k)$ coordinates system. 
     360 
     361%% ================================================================================================= 
    382362\subsection{Continuous model equations} 
    383 \label{subsec:PE_zco_Eq} 
     363\label{subsec:MB_zco_Eq} 
    384364 
    385365In order to express the Primitive Equations in tensorial formalism, 
    386 it is necessary to compute the horizontal component of the non-linear and viscous terms of the equation using 
    387 \autoref{eq:PE_grad}) to \autoref{eq:PE_lap_vector}. 
    388 Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, the velocity in the $(i,j,k)$ coordinates system and 
    389 define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, by: 
     366it is necessary to compute the horizontal component of the non-linear and viscous terms of 
     367the equation using \autoref{eq:MB_grad}) to \autoref{eq:MB_lap_vector}. 
     368Let us set $\vect U = (u,v,w) = \vect U_h + w \; \vect k $, 
     369the velocity in the $(i,j,k)$ coordinates system, 
     370and define the relative vorticity $\zeta$ and the divergence of the horizontal velocity field $\chi$, 
     371by: 
    390372\begin{gather} 
    391   \label{eq:PE_curl_Uh} 
     373  \label{eq:MB_curl_Uh} 
    392374  \zeta = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, v)]{i} - \pd[(e_1 \, u)]{j} \rt] \\ 
    393   \label{eq:PE_div_Uh} 
     375  \label{eq:MB_div_Uh} 
    394376  \chi  = \frac{1}{e_1 e_2} \lt[ \pd[(e_2 \, u)]{i} + \pd[(e_1 \, v)]{j} \rt] 
    395377\end{gather} 
     
    397379Using again the fact that the horizontal scale factors $e_1$ and $e_2$ are independent of $k$ and that 
    398380$e_3$  is a function of the single variable $k$, 
    399 $NLT$ the nonlinear term of \autoref{eq:PE_dyn} can be transformed as follows: 
    400 \begin{alignat*}{2} 
    401   &NLT &=   &\lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ 
    402   &    &=   &\lt( 
    403     \begin{array}{*{20}c} 
    404                 \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v   \\ 
    405                 \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w 
    406     \end{array} 
    407                                                                                              \rt) 
    408           + \frac{1}{2} \lt( 
    409     \begin{array}{*{20}c} 
    410                              \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ 
    411                              \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} 
    412     \end{array} 
    413                                                                      \rt) \\ 
    414   &    &=   &\lt( 
    415     \begin{array}{*{20}c} 
    416                   -\zeta \; v \\ 
    417                    \zeta \; u 
    418     \end{array} 
    419                               \rt) 
    420           + \frac{1}{2} \lt( 
    421     \begin{array}{*{20}c} 
    422                              \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ 
    423                              \frac{1}{e_2} \pd[(u^2 + v^2)]{j} 
    424     \end{array} 
    425                                                                \rt) \\ 
    426   &    &  &+ \frac{1}{e_3} \lt( 
    427     \begin{array}{*{20}c} 
    428                                 w \; \pd[u]{k} \\ 
    429                                 w \; \pd[v]{k} 
    430     \end{array} 
    431                                                \rt) 
    432            - \lt( 
    433     \begin{array}{*{20}c} 
    434                   \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\ 
    435                   \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j} 
    436     \end{array} 
    437                                                                         \rt) 
    438 \end{alignat*} 
    439 The last term of the right hand side is obviously zero, and thus the nonlinear term of 
    440 \autoref{eq:PE_dyn} is written in the $(i,j,k)$ coordinate system: 
     381$NLT$ the nonlinear term of \autoref{eq:MB_PE_dyn} can be transformed as follows: 
     382\begin{align*} 
     383  NLT &= \lt[ (\nabla \times {\vect U}) \times {\vect U} + \frac{1}{2} \nabla \lt( {\vect U}^2 \rt) \rt]_h \\ 
     384      &= \lt( 
     385        \begin{array}{*{20}c} 
     386          \lt[ \frac{1}{e_3} \pd[u]{k} - \frac{1}{e_1} \pd[w]{i} \rt] w - \zeta \; v   \\ 
     387          \zeta \; u - \lt[ \frac{1}{e_2} \pd[w]{j} - \frac{1}{e_3} \pd[v]{k} \rt] \ w 
     388        \end{array} 
     389  \rt) 
     390  + \frac{1}{2} \lt( 
     391  \begin{array}{*{20}c} 
     392    \frac{1}{e_1} \pd[(u^2 + v^2 + w^2)]{i} \\ 
     393    \frac{1}{e_2} \pd[(u^2 + v^2 + w^2)]{j} 
     394  \end{array} 
     395  \rt) \\ 
     396      &= \lt( 
     397        \begin{array}{*{20}c} 
     398          -\zeta \; v \\ 
     399          \zeta \; u 
     400        \end{array} 
     401  \rt) 
     402  + \frac{1}{2} \lt( 
     403  \begin{array}{*{20}c} 
     404    \frac{1}{e_1} \pd[(u^2 + v^2)]{i} \\ 
     405    \frac{1}{e_2} \pd[(u^2 + v^2)]{j} 
     406  \end{array} 
     407  \rt) 
     408  + \frac{1}{e_3} \lt( 
     409  \begin{array}{*{20}c} 
     410    w \; \pd[u]{k} \\ 
     411    w \; \pd[v]{k} 
     412  \end{array} 
     413  \rt) 
     414  - \lt( 
     415  \begin{array}{*{20}c} 
     416    \frac{w}{e_1} \pd[w]{i} - \frac{1}{2 e_1} \pd[w^2]{i} \\ 
     417    \frac{w}{e_2} \pd[w]{j} - \frac{1}{2 e_2} \pd[w^2]{j} 
     418  \end{array} 
     419  \rt) 
     420\end{align*} 
     421The last term of the right hand side is obviously zero, 
     422and thus the \textbf{N}on\textbf{L}inear \textbf{T}erm ($NLT$) of \autoref{eq:MB_PE_dyn} is written in 
     423the $(i,j,k)$ coordinate system: 
    441424\begin{equation} 
    442   \label{eq:PE_vector_form} 
    443   NLT =   \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) 
    444         + \frac{1}{e_3} w \pd[\vect U_h]{k} 
     425  \label{eq:MB_vector_form} 
     426  NLT = \zeta \; \vect k \times \vect U_h + \frac{1}{2} \nabla_h \lt( \vect U_h^2 \rt) + \frac{1}{e_3} w \pd[\vect U_h]{k} 
    445427\end{equation} 
    446428 
     
    448430For some purposes, it can be advantageous to write this term in the so-called flux form, 
    449431\ie\ to write it as the divergence of fluxes. 
    450 For example, the first component of \autoref{eq:PE_vector_form} (the $i$-component) is transformed as follows: 
    451 \begin{alignat*}{2} 
     432For example, 
     433the first component of \autoref{eq:MB_vector_form} (the $i$-component) is transformed as follows: 
     434\begin{alignat*}{3} 
    452435  &NLT_i &= &- \zeta \; v + \frac{1}{2 \; e_1} \pd[ (u^2 + v^2) ]{i} + \frac{1}{e_3} w \ \pd[u]{k} \\ 
    453   &      &=  &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) 
    454             + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) \\ 
    455   &      & &+ \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ 
    456   &      &=  &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) 
    457                                      + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} -         e_1 \, u \pd[v]{j} \rt) \rt. \\ 
    458   &      &                       &\lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt) 
    459                                      + e_2 v \pd[v]{i}                                                         \rt] \\ 
     436  &      &= &\frac{1}{e_1 \; e_2} \lt( -v \pd[(e_2 \, v)]{i} + v \pd[(e_1 \, u)]{j} \rt) + \frac{1}{e_1 e_2} \lt( e_2 \; u \pd[u]{i} + e_2 \; v \pd[v]{i} \rt) + \frac{1}{e_3} \lt( w \; \pd[u]{k} \rt) \\ 
     437  &      &= &\frac{1}{e_1 \; e_2} \lt[ - \lt( v^2 \pd[e_2]{i} + e_2 \, v \pd[v]{i} \rt) + \lt( \pd[ \lt( e_1 \, u \, v \rt)]{j} - e_1 \, u \pd[v]{j} \rt) \rt. \lt. + \lt( \pd[ \lt( e_2 \, u \, u \rt)]{i} - u \pd[ \lt( e_2 u \rt)]{i} \rt) + e_2 v \pd[v]{i} \rt] \\ 
    460438  &      & &+ \frac{1}{e_3} \lt( \pd[(w \, u)]{k} - u \pd[w]{k} \rt) \\ 
    461   &      &=  &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) 
    462             + \frac{1}{e_3} \pd[(w \, u)]{k} \\ 
    463   &      & &+ \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt) 
    464                                   - u \pd[(e_2 u)]{i}                              \rt] 
    465             - \frac{1}{e_3} \pd[w]{k} u \\ 
     439  &      &= &\frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, u \, v)]{j} \rt) + \frac{1}{e_3} \pd[(w \, u)]{k} + \frac{1}{e_1 e_2} \lt[ - u \lt( \pd[(e_1 v)]{j} - v \, \pd[e_1]{j} \rt) - u \pd[(e_2 u)]{i} \rt] - \frac{1}{e_3} \pd[w]{k} u \\ 
    466440  &      & &+ \frac{1}{e_1 e_2} \lt( - v^2 \pd[e_2]{i} \rt) \\ 
    467   &      &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u 
    468             + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 
    469   \intertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} 
     441  &      &= &\nabla \cdot (\vect U \, u) - (\nabla \cdot \vect U) \ u + \frac{1}{e_1 e_2} \lt( -v^2 \pd[e_2]{i} + u v \, \pd[e_1]{j} \rt) \\ 
     442  \shortintertext{as $\nabla \cdot {\vect U} \; = 0$ (incompressibility) it becomes:} 
    470443  &      &= &\, \nabla \cdot (\vect U \, u) + \frac{1}{e_1 e_2} \lt( v \; \pd[e_2]{i} - u \; \pd[e_1]{j} \rt) (-v) 
    471444\end{alignat*} 
     
    473446The flux form of the momentum advection term is therefore given by: 
    474447\begin{equation} 
    475   \label{eq:PE_flux_form} 
    476   NLT =   \nabla \cdot \lt( 
    477     \begin{array}{*{20}c} 
    478                             \vect U \, u \\ 
    479                             \vect U \, v 
    480     \end{array} 
    481                                          \rt) 
    482         + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h 
     448  \label{eq:MB_flux_form} 
     449  NLT = \nabla \cdot \lt( 
     450  \begin{array}{*{20}c} 
     451    \vect U \, u \\ 
     452    \vect U \, v 
     453  \end{array} 
     454  \rt) 
     455  + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \vect k \times \vect U_h 
    483456\end{equation} 
    484457 
    485458The flux form has two terms, 
    486 the first one is expressed as the divergence of momentum fluxes (hence the flux form name given to this formulation) 
    487 and the second one is due to the curvilinear nature of the coordinate system used. 
    488 The latter is called the \textit{metric} term and can be viewed as a modification of the Coriolis parameter: 
     459the first one is expressed as the divergence of momentum fluxes 
     460(hence the flux form name given to this formulation) and 
     461the second one is due to the curvilinear nature of the coordinate system used. 
     462The latter is called the \textit{metric} term and can be viewed as 
     463a modification of the Coriolis parameter: 
    489464\[ 
    490   % \label{eq:PE_cor+metric} 
     465  % \label{eq:MB_cor+metric} 
    491466  f \to f + \frac{1}{e_1 e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) 
    492467\] 
     
    499474the following tensorial formalism: 
    500475 
    501 \begin{itemize} 
    502 \item 
    503   \textbf{Vector invariant form of the momentum equations}: 
     476\begin{description} 
     477\item [Vector invariant form of the momentum equations] 
    504478  \begin{equation} 
    505     \label{eq:PE_dyn_vect} 
    506     \begin{split} 
    507     % \label{eq:PE_dyn_vect_u} 
    508       \pd[u]{t} = &+ (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) 
    509                    - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 
    510                   &+ D_u^{\vect U} + F_u^{\vect U} \\ 
    511       \pd[v]{t} = &- (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2) 
    512                    - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) \\ 
    513                   &+ D_v^{\vect U} + F_v^{\vect U} 
    514     \end{split} 
     479    \label{eq:MB_dyn_vect} 
     480    \begin{gathered} 
     481    % \label{eq:MB_dyn_vect_u} 
     482      \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} w \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} \\ 
     483      \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 e_2} \pd[]{j} (u^2 + v^2) - \frac{1}{e_3} w \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} 
     484    \end{gathered} 
    515485  \end{equation} 
    516 \item 
    517   \textbf{flux form of the momentum equations}: 
    518   % \label{eq:PE_dyn_flux} 
    519   \begin{multline*} 
    520     % \label{eq:PE_dyn_flux_u} 
    521     \pd[u]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\ 
    522                 - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ 
    523                 - \frac{1}{e_3} \pd[(w \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
    524                 + D_u^{\vect U} + F_u^{\vect U} 
    525   \end{multline*} 
    526   \begin{multline*} 
    527     % \label{eq:PE_dyn_flux_v} 
    528     \pd[v]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\ 
    529                 - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 
    530                 - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
    531                 + D_v^{\vect U} + F_v^{\vect U} 
    532   \end{multline*} 
    533   where $\zeta$, the relative vorticity, is given by \autoref{eq:PE_curl_Uh} and $p_s$, the surface pressure, 
    534   is given by: 
     486\item [Flux form of the momentum equations] 
     487  % \label{eq:MB_dyn_flux} 
     488  \begin{alignat*}{2} 
     489    % \label{eq:MB_dyn_flux_u} 
     490    \pd[u]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(w \, u)]{k} \\ 
     491    &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} 
     492  \end{alignat*} 
     493  \begin{alignat*}{2} 
     494    % \label{eq:MB_dyn_flux_v} 
     495    \pd[v]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(w \, v)]{k} \\ 
     496    &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} 
     497  \end{alignat*} 
     498  where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and 
     499  $p_s$, the surface pressure, is given by: 
    535500  \[ 
    536   % \label{eq:PE_spg} 
     501  % \label{eq:MB_spg} 
    537502    p_s = \rho \,g \, \eta 
    538503  \] 
    539   and $\eta$ is the solution of \autoref{eq:PE_ssh}. 
     504  and $\eta$ is the solution of \autoref{eq:MB_ssh}. 
    540505 
    541506  The vertical velocity and the hydrostatic pressure are diagnosed from the following equations: 
    542507  \[ 
    543   % \label{eq:w_diag} 
     508  % \label{eq:MB_w_diag} 
    544509    \pd[w]{k} = - \chi \; e_3 \qquad 
    545   % \label{eq:hp_diag} 
     510  % \label{eq:MB_hp_diag} 
    546511    \pd[p_h]{k} = - \rho \; g \; e_3 
    547512  \] 
    548   where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:PE_div_Uh}. 
    549  
    550 \item 
    551   \textbf{tracer equations}: 
    552   \begin{equation} 
    553   \begin{split} 
    554     \pd[T]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] 
    555                 - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 
    556     \pd[S]{t} = & - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] 
    557                 - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ 
    558     \rho = & \rho \big( T,S,z(k) \big) 
    559   \end{split} 
    560   \end{equation} 
    561 \end{itemize} 
    562  
    563 The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on the subgrid scale parameterisation used. 
    564 It will be defined in \autoref{eq:PE_zdf}. 
     513  where the divergence of the horizontal velocity, $\chi$ is given by \autoref{eq:MB_div_Uh}. 
     514\item [Tracer equations] 
     515  \begin{gather*} 
     516    \pd[T]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 T \, u)]{i} + \pd[(e_1 T \, v)]{j} \rt] - \frac{1}{e_3} \pd[(T \, w)]{k} + D^T + F^T \\ 
     517    \pd[S]{t} = - \frac{1}{e_1 e_2} \lt[ \pd[(e_2 S \, u)]{i} + \pd[(e_1 S \, v)]{j} \rt] - \frac{1}{e_3} \pd[(S \, w)]{k} + D^S + F^S \\ 
     518    \rho = \rho \big( T,S,z(k) \big) 
     519  \end{gather*} 
     520\end{description} 
     521 
     522The expression of $\vect D^{U}$, $D^{S}$ and $D^{T}$ depends on 
     523the subgrid scale parameterisation used. 
     524It will be defined in \autoref{eq:MB_zdf}. 
    565525The nature and formulation of $\vect F^{\vect U}$, $F^T$ and $F^S$, the surface forcing terms, 
    566526are discussed in \autoref{chap:SBC}. 
    567527 
    568 \newpage 
    569  
    570 % ================================================================ 
    571 % Curvilinear generalised vertical coordinate System 
    572 % ================================================================ 
     528%% ================================================================================================= 
    573529\section{Curvilinear generalised vertical coordinate system} 
    574 \label{sec:PE_gco} 
     530\label{sec:MB_gco} 
    575531 
    576532The ocean domain presents a huge diversity of situation in the vertical. 
     
    579535varying from more than 6,000 meters in abyssal trenches to zero at the coast. 
    580536Last but not least, the ocean stratification exerts a strong barrier to vertical motions and mixing. 
    581 Therefore, in order to represent the ocean with respect to 
    582 the first point a space and time dependent vertical coordinate that follows the variation of the sea surface height 
     537Therefore, in order to represent the ocean with respect to the first point 
     538a space and time dependent vertical coordinate that follows the variation of the sea surface height 
    583539\eg\ an \zstar-coordinate; 
    584 for the second point, a space variation to fit the change of bottom topography 
     540for the second point, 
     541a space variation to fit the change of bottom topography 
    585542\eg\ a terrain-following or $\sigma$-coordinate; 
    586 and for the third point, one will be tempted to use a space and time dependent coordinate that 
    587 follows the isopycnal surfaces, \eg\ an isopycnic coordinate. 
    588  
    589 In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, as in 
    590 HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and $\sigma$ at 
    591 the ocean bottom) \citep{chassignet.smith.ea_JPO03} or 
    592 OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and $\sigma$-coordinate elsewhere) 
    593 \citep{madec.delecluse.ea_JPO96} among others. 
     543and for the third point, 
     544one will be tempted to use a space and time dependent coordinate that follows the isopycnal surfaces, 
     545\eg\ an isopycnic coordinate. 
     546 
     547In order to satisfy two or more constraints one can even be tempted to mixed these coordinate systems, 
     548as in HYCOM (mixture of $z$-coordinate at the surface, isopycnic coordinate in the ocean interior and 
     549$\sigma$ at the ocean bottom) \citep{chassignet.smith.ea_JPO03} or 
     550OPA (mixture of $z$-coordinate in vicinity the surface and steep topography areas and 
     551$\sigma$-coordinate elsewhere) \citep{madec.delecluse.ea_JPO96} among others. 
    594552 
    595553In fact one is totally free to choose any space and time vertical coordinate by 
    596554introducing an arbitrary vertical coordinate : 
    597555\begin{equation} 
    598   \label{eq:PE_s} 
     556  \label{eq:MB_s} 
    599557  s = s(i,j,k,t) 
    600558\end{equation} 
    601 with the restriction that the above equation gives a single-valued monotonic relationship between $s$ and $k$, 
    602 when $i$, $j$ and $t$ are held fixed. 
    603 \autoref{eq:PE_s} is a transformation from the $(i,j,k,t)$ coordinate system with independent variables into 
     559with the restriction that the above equation gives a single-valued monotonic relationship between 
     560$s$ and $k$, when $i$, $j$ and $t$ are held fixed. 
     561\autoref{eq:MB_s} is a transformation from 
     562the $(i,j,k,t)$ coordinate system with independent variables into 
    604563the $(i,j,s,t)$ generalised coordinate system with $s$ depending on the other three variables through 
    605 \autoref{eq:PE_s}. 
     564\autoref{eq:MB_s}. 
    606565This so-called \textit{generalised vertical coordinate} \citep{kasahara_MWR74} is in fact 
    607 an Arbitrary Lagrangian--Eulerian (ALE) coordinate. 
    608 Indeed, one has a great deal of freedom in the choice of expression for $s$. The choice determines 
    609 which part of the vertical velocity (defined from a fixed referential) will cross the levels (Eulerian part) and 
    610 which part will be used to move them (Lagrangian part). 
    611 The coordinate is also sometime referenced as an adaptive coordinate \citep{hofmeister.burchard.ea_OM10}, 
    612 since the coordinate system is adapted in the course of the simulation. 
     566an \textbf{A}rbitrary \textbf{L}agrangian--\textbf{E}ulerian (ALE) coordinate. 
     567Indeed, one has a great deal of freedom in the choice of expression for $s$. 
     568The choice determines which part of the vertical velocity (defined from a fixed referential) 
     569will cross the levels (Eulerian part) and which part will be used to move them (Lagrangian part). 
     570The coordinate is also sometimes referenced as an adaptive coordinate 
     571\citep{hofmeister.burchard.ea_OM10}, since the coordinate system is adapted in 
     572the course of the simulation. 
    613573Its most often used implementation is via an ALE algorithm, 
    614574in which a pure lagrangian step is followed by regridding and remapping steps, 
     
    627587The horizontal velocity in ocean models measures motions in the horizontal plane, 
    628588perpendicular to the local gravitational field. 
    629 That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, be it geopotential, 
    630 isopycnal, pressure, or terrain following. 
     589That is, horizontal velocity is mathematically the same regardless of the vertical coordinate, 
     590be it geopotential, isopycnal, pressure, or terrain following. 
    631591The key motivation for maintaining the same horizontal velocity component is that 
    632592the hydrostatic and geostrophic balances are dominant in the large-scale ocean. 
    633 Use of an alternative quasi -horizontal velocity, for example one oriented parallel to the generalized surface, 
     593Use of an alternative quasi -horizontal velocity, 
     594for example one oriented parallel to the generalized surface, 
    634595would lead to unacceptable numerical errors. 
    635596Correspondingly, the vertical direction is anti -parallel to the gravitational force in 
     
    638599the surface of a constant generalized vertical coordinate. 
    639600 
    640 It is the method used to measure transport across the generalized vertical coordinate surfaces which differs between 
    641 the vertical coordinate choices. 
    642 That is, computation of the dia-surface velocity component represents the fundamental distinction between 
     601It is the method used to measure transport across the generalized vertical coordinate surfaces which 
     602differs between the vertical coordinate choices. 
     603That is, 
     604computation of the dia-surface velocity component represents the fundamental distinction between 
    643605the various coordinates. 
    644 In some models, such as geopotential, pressure, and terrain following, this transport is typically diagnosed from 
    645 volume or mass conservation. 
    646 In other models, such as isopycnal layered models, this transport is prescribed based on assumptions about 
    647 the physical processes producing a flux across the layer interfaces. 
     606In some models, such as geopotential, pressure, and terrain following, 
     607this transport is typically diagnosed from volume or mass conservation. 
     608In other models, such as isopycnal layered models, 
     609this transport is prescribed based on assumptions about the physical processes producing 
     610a flux across the layer interfaces. 
    648611 
    649612In this section we first establish the PE in the generalised vertical $s$-coordinate, 
     
    651614%} 
    652615 
    653 % ------------------------------------------------------------------------------------------------------------- 
    654 % The s-coordinate Formulation 
    655 % ------------------------------------------------------------------------------------------------------------- 
     616%% ================================================================================================= 
    656617\subsection{\textit{S}-coordinate formulation} 
    657618 
    658 Starting from the set of equations established in \autoref{sec:PE_zco} for the special case $k = z$ and 
    659 thus $e_3 = 1$, we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, 
     619Starting from the set of equations established in \autoref{sec:MB_zco} for 
     620the special case $k = z$ and thus $e_3 = 1$, 
     621we introduce an arbitrary vertical coordinate $s = s(i,j,k,t)$, 
    660622which includes $z$-, \zstar- and $\sigma$-coordinates as special cases 
    661623($s = z$, $s = \zstar$, and $s = \sigma = z / H$ or $ = z / \lt( H + \eta \rt)$, resp.). 
    662 A formal derivation of the transformed equations is given in \autoref{apdx:A}. 
    663 Let us define the vertical scale factor by $e_3 = \partial_s z$  ($e_3$ is now a function of $(i,j,k,t)$ ), 
     624A formal derivation of the transformed equations is given in \autoref{apdx:SCOORD}. 
     625Let us define the vertical scale factor by $e_3 = \partial_s z$ 
     626($e_3$ is now a function of $(i,j,k,t)$ ), 
    664627and the slopes in the $(i,j)$ directions between $s$- and $z$-surfaces by: 
    665628\begin{equation} 
    666   \label{eq:PE_sco_slope} 
     629  \label{eq:MB_sco_slope} 
    667630  \sigma_1 = \frac{1}{e_1} \; \lt. \pd[z]{i} \rt|_s \quad \text{and} \quad 
    668631  \sigma_2 = \frac{1}{e_2} \; \lt. \pd[z]{j} \rt|_s 
    669632\end{equation} 
    670 We also introduce $\omega$, a dia-surface velocity component, defined as the velocity 
    671 relative to the moving $s$-surfaces and normal to them: 
     633We also introduce $\omega$, a dia-surface velocity component, 
     634defined as the velocity relative to the moving $s$-surfaces and normal to them: 
    672635\[ 
    673   % \label{eq:PE_sco_w} 
     636  % \label{eq:MB_sco_w} 
    674637  \omega = w -  \, \lt. \pd[z]{t} \rt|_s - \sigma_1 \, u - \sigma_2 \, v 
    675638\] 
    676639 
    677 The equations solved by the ocean model \autoref{eq:PE} in $s$-coordinate can be written as follows 
    678 (see \autoref{sec:A_momentum}): 
    679  
    680 \begin{itemize} 
    681 \item \textbf{Vector invariant form of the momentum equation}: 
    682   \begin{multline*} 
    683   % \label{eq:PE_sco_u_vector} 
    684     \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} \\ 
    685                 - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 
    686                 + D_u^{\vect U} + F_u^{\vect U} 
    687   \end{multline*} 
    688   \begin{multline*} 
    689   % \label{eq:PE_sco_v_vector} 
    690     \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j}(u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} \\ 
    691                 - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_2 
    692                 + D_v^{\vect U} + F_v^{\vect U} 
    693   \end{multline*} 
    694 \item \textbf{Flux form of the momentum equation}: 
    695   \begin{multline*} 
    696   % \label{eq:PE_sco_u_flux} 
    697     \frac{1}{e_3} \pd[(e_3 \, u)]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v \\ 
    698                                        - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) \\ 
    699                                        - \frac{1}{e_3} \pd[(\omega \, u)]{k} 
    700                                        - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
    701                                        - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 
    702   \end{multline*} 
    703   \begin{multline*} 
    704   % \label{eq:PE_sco_v_flux} 
    705     \frac{1}{e_3} \pd[(e_3 \, v)]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u \\ 
    706                                        - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) \\ 
    707                                        - \frac{1}{e_3} \pd[(\omega \, v)]{k} 
    708                                        - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) 
    709                                        - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
    710   \end{multline*} 
     640The equations solved by the ocean model \autoref{eq:MB_PE} in $s$-coordinate can be written as follows 
     641(see \autoref{sec:SCOORD_momentum}): 
     642 
     643\begin{description} 
     644\item [Vector invariant form of the momentum equation] 
     645  \begin{gather*} 
     646    % \label{eq:MB_sco_u_vector} 
     647    \pd[u]{t} = + (\zeta + f) \, v - \frac{1}{2 \, e_1} \pd[]{i} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[u]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} \\ 
     648    % \label{eq:MB_sco_v_vector} 
     649    \pd[v]{t} = - (\zeta + f) \, u - \frac{1}{2 \, e_2} \pd[]{j} (u^2 + v^2) - \frac{1}{e_3} \omega \pd[v]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
     650  \end{gather*} 
     651\item [Flux form of the momentum equation] 
     652  \begin{alignat*}{2} 
     653    % \label{eq:MB_sco_u_flux} 
     654    \frac{1}{e_3} \pd[(e_3 \, u)]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, u)]{k} \\ 
     655    &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 
     656  \end{alignat*} 
     657  \begin{alignat*}{2} 
     658  % \label{eq:MB_sco_v_flux} 
     659    \frac{1}{e_3} \pd[(e_3 \, v)]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, v)]{k} \\ 
     660    &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
     661  \end{alignat*} 
    711662  where the relative vorticity, $\zeta$, the surface pressure gradient, 
    712663  and the hydrostatic pressure have the same expressions as in $z$-coordinates although 
    713664  they do not represent exactly the same quantities. 
    714   $\omega$ is provided by the continuity equation (see \autoref{apdx:A}): 
     665  $\omega$ is provided by the continuity equation (see \autoref{apdx:SCOORD}): 
    715666  \[ 
    716   % \label{eq:PE_sco_continuity} 
     667  % \label{eq:MB_sco_continuity} 
    717668    \pd[e_3]{t} + e_3 \; \chi + \pd[\omega]{s} = 0 \quad \text{with} \quad 
    718669    \chi = \frac{1}{e_1 e_2 e_3} \lt( \pd[(e_2 e_3 \, u)]{i} + \pd[(e_1 e_3 \, v)]{j} \rt) 
    719670  \] 
    720 \item \textit{tracer equations}: 
    721   \begin{multline*} 
    722   % \label{eq:PE_sco_t} 
    723     \frac{1}{e_3} \pd[(e_3 \, T)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, T)]{i} 
    724                                                                     + \pd[(e_1 e_3 \, v \, T)]{j} \rt) \\ 
    725                                        - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S 
    726   \end{multline*} 
    727   \begin{multline} 
    728   % \label{eq:PE_sco_s} 
    729     \frac{1}{e_3} \pd[(e_3 \, S)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, S)]{i} 
    730                                                                     + \pd[(e_1 e_3 \, v \, S)]{j} \rt) \\ 
    731                                        - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S 
    732   \end{multline} 
    733 \end{itemize} 
     671\item [Tracer equations] 
     672  \begin{gather*} 
     673    % \label{eq:MB_sco_t} 
     674    \frac{1}{e_3} \pd[(e_3 \, T)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, T)]{i} + \pd[(e_1 e_3 \, v \, T)]{j} \rt) - \frac{1}{e_3} \pd[(T \, \omega)]{k} + D^T + F^S \\ 
     675    % \label{eq:MB_sco_s} 
     676    \frac{1}{e_3} \pd[(e_3 \, S)]{t} = - \frac{1}{e_1 e_2 e_3} \lt(   \pd[(e_2 e_3 \, u \, S)]{i} + \pd[(e_1 e_3 \, v \, S)]{j} \rt) - \frac{1}{e_3} \pd[(S \, \omega)]{k} + D^S + F^S 
     677  \end{gather*} 
     678\end{description} 
    734679The equation of state has the same expression as in $z$-coordinate, 
    735680and similar expressions are used for mixing and forcing terms. 
     
    741686} 
    742687 
    743 % ------------------------------------------------------------------------------------------------------------- 
    744 % Curvilinear \zstar-coordinate System 
    745 % ------------------------------------------------------------------------------------------------------------- 
     688%% ================================================================================================= 
    746689\subsection{Curvilinear \zstar-coordinate system} 
    747 \label{subsec:PE_zco_star} 
    748  
    749 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    750 \begin{figure}[!b] 
    751   \begin{center} 
    752     \includegraphics[width=\textwidth]{Fig_z_zstar} 
    753     \caption{ 
    754       \protect\label{fig:z_zstar} 
    755       (a) $z$-coordinate in linear free-surface case ; 
    756       (b) $z$-coordinate in non-linear free surface case ; 
    757       (c) re-scaled height coordinate 
     690\label{subsec:MB_zco_star} 
     691 
     692\begin{figure} 
     693  \centering 
     694  \includegraphics[width=0.66\textwidth]{Fig_z_zstar} 
     695  \caption[Curvilinear z-coordinate systems (\{non-\}linear free-surface cases and re-scaled \zstar)]{ 
     696    \begin{enumerate*}[label=(\textit{\alph*})] 
     697    \item $z$-coordinate in linear free-surface case; 
     698    \item $z$-coordinate in non-linear free surface case; 
     699    \item re-scaled height coordinate 
    758700      (become popular as the \zstar-coordinate \citep{adcroft.campin_OM04}). 
    759     } 
    760   \end{center} 
     701    \end{enumerate*} 
     702  } 
     703  \label{fig:MB_z_zstar} 
    761704\end{figure} 
    762 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    763  
    764 In this case, the free surface equation is nonlinear, and the variations of volume are fully taken into account. 
    765 These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on the \NEMO\ web site. 
    766  
    767 The \zstar coordinate approach is an unapproximated, non-linear free surface implementation which allows one to 
    768 deal with large amplitude free-surface variations relative to the vertical resolution \citep{adcroft.campin_OM04}. 
    769 In the \zstar formulation, 
    770 the variation of the column thickness due to sea-surface undulations is not concentrated in the surface level, 
    771 as in the $z$-coordinate formulation, but is equally distributed over the full water column. 
     705 
     706In this case, the free surface equation is nonlinear, 
     707and the variations of volume are fully taken into account. 
     708These coordinates systems is presented in a report \citep{levier.treguier.ea_rpt07} available on 
     709the \NEMO\ web site. 
     710 
     711The \zstar\ coordinate approach is an unapproximated, non-linear free surface implementation which 
     712allows one to deal with large amplitude free-surface variations relative to the vertical resolution 
     713\citep{adcroft.campin_OM04}. 
     714In the \zstar\ formulation, the variation of the column thickness due to sea-surface undulations is 
     715not concentrated in the surface level, as in the $z$-coordinate formulation, 
     716but is equally distributed over the full water column. 
    772717Thus vertical levels naturally follow sea-surface variations, with a linear attenuation with depth, 
    773 as illustrated by \autoref{fig:z_zstar}. 
    774 Note that with a flat bottom, such as in \autoref{fig:z_zstar}, the bottom-following $z$ coordinate and \zstar are equivalent. 
     718as illustrated by \autoref{fig:MB_z_zstar}. 
     719Note that with a flat bottom, such as in \autoref{fig:MB_z_zstar}, 
     720the bottom-following $z$ coordinate and \zstar\ are equivalent. 
    775721The definition and modified oceanic equations for the rescaled vertical coordinate \zstar, 
    776722including the treatment of fresh-water flux at the surface, are detailed in Adcroft and Campin (2004). 
    777723The major points are summarized here. 
    778 The position (\zstar) and vertical discretization (\zstar) are expressed as: 
     724The position (\zstar) and vertical discretization ($\delta \zstar$) are expressed as: 
    779725\[ 
    780   % \label{eq:PE_z-star} 
    781   H + \zstar = (H + z)  / r \quad \text{and}  \quad \delta \zstar 
    782               = \delta z / r \quad \text{with} \quad r 
    783               = \frac{H + \eta}{H} . 
     726  % \label{eq:MB_z-star} 
     727  H + \zstar = (H + z)  / r \quad \text{and}  \quad \delta \zstar = \delta z / r \quad \text{with} \quad r = \frac{H + \eta}{H} 
    784728\] 
    785729Simple re-organisation of the above expressions gives 
    786730\[ 
    787   % \label{eq:PE_zstar_2} 
    788   \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) . 
     731  % \label{eq:MB_zstar_2} 
     732  \zstar = H \lt( \frac{z - \eta}{H + \eta} \rt) 
    789733\] 
    790734Since the vertical displacement of the free surface is incorporated in the vertical coordinate \zstar, 
    791 the upper and lower boundaries are at fixed  \zstar position, 
     735the upper and lower boundaries are at fixed \zstar\ position, 
    792736$\zstar = 0$ and $\zstar = -H$ respectively. 
    793737Also the divergence of the flow field is no longer zero as shown by the continuity equation: 
    794738\[ 
    795   \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 . 
     739  \pd[r]{t} = \nabla_{\zstar} \cdot \lt( r \; \vect U_h \rt) + \pd[r \; w^*]{\zstar} = 0 
    796740\] 
    797 This \zstar coordinate is closely related to the "eta" coordinate used in many atmospheric models 
    798 (see Black (1994) for a review of eta coordinate atmospheric models). 
     741This \zstar\ coordinate is closely related to the $\eta$ coordinate used in many atmospheric models 
     742(see Black (1994) for a review of $\eta$ coordinate atmospheric models). 
    799743It was originally used in ocean models by Stacey et al. (1995) for studies of tides next to shelves, 
    800744and it has been recently promoted by Adcroft and Campin (2004) for global climate modelling. 
    801745 
    802 The surfaces of constant \zstar are quasi -horizontal. 
    803 Indeed, the \zstar coordinate reduces to $z$ when $\eta$ is zero. 
     746The surfaces of constant \zstar\ are quasi-horizontal. 
     747Indeed, the \zstar\ coordinate reduces to $z$ when $\eta$ is zero. 
    804748In general, when noting the large differences between 
    805749undulations of the bottom topography versus undulations in the surface height, 
    806 it is clear that surfaces constant \zstar are very similar to the depth surfaces. 
     750it is clear that surfaces constant \zstar\ are very similar to the depth surfaces. 
    807751These properties greatly reduce difficulties of computing the horizontal pressure gradient relative to 
    808 terrain following sigma models discussed in \autoref{subsec:PE_sco}. 
     752terrain following sigma models discussed in \autoref{subsec:MB_sco}. 
    809753Additionally, since $\zstar = z$ when $\eta = 0$, 
    810 no flow is spontaneously generated in an unforced ocean starting from rest, regardless the bottom topography. 
    811 This behaviour is in contrast to the case with "s"-models, where pressure gradient errors in the presence of 
    812 nontrivial topographic variations can generate nontrivial spontaneous flow from a resting state, 
     754no flow is spontaneously generated in an unforced ocean starting from rest, 
     755regardless the bottom topography. 
     756This behaviour is in contrast to the case with ``s''-models, 
     757where pressure gradient errors in the presence of nontrivial topographic variations can generate 
     758nontrivial spontaneous flow from a resting state, 
    813759depending on the sophistication of the pressure gradient solver. 
    814 The quasi -horizontal nature of the coordinate surfaces also facilitates the implementation of 
    815 neutral physics parameterizations in \zstar  models using the same techniques as in $z$-models 
     760The quasi-horizontal nature of the coordinate surfaces also facilitates the implementation of 
     761neutral physics parameterizations in \zstar\ models using the same techniques as in $z$-models 
    816762(see Chapters 13-16 of \cite{griffies_bk04}) for a discussion of neutral physics in $z$-models, 
    817763as well as \autoref{sec:LDF_slp} in this document for treatment in \NEMO). 
    818764 
    819 The range over which \zstar  varies is time independent $-H \leq \zstar \leq 0$. 
     765The range over which \zstar\ varies is time independent $-H \leq \zstar \leq 0$. 
    820766Hence, all cells remain nonvanishing, so long as the surface height maintains $\eta > -H$. 
    821 This is a minor constraint relative to that encountered on the surface height when using $s = z$ or $s = z - \eta$. 
    822  
    823 Because \zstar  has a time independent range, all grid cells have static increments ds, 
     767This is a minor constraint relative to that encountered on the surface height when 
     768using $s = z$ or $s = z - \eta$. 
     769 
     770Because \zstar\ has a time independent range, all grid cells have static increments ds, 
    824771and the sum of the vertical increments yields the time independent ocean depth. %k ds = H. 
    825 The \zstar coordinate is therefore invisible to undulations of the free surface, 
     772The \zstar\ coordinate is therefore invisible to undulations of the free surface, 
    826773since it moves along with the free surface. 
    827 This property means that no spurious vertical transport is induced across surfaces of constant \zstar  by 
    828 the motion of external gravity waves. 
     774This property means that no spurious vertical transport is induced across 
     775surfaces of constant \zstar\  by the motion of external gravity waves. 
    829776Such spurious transport can be a problem in z-models, especially those with tidal forcing. 
    830 Quite generally, the time independent range for the \zstar  coordinate is a very convenient property that 
    831 allows for a nearly arbitrary vertical resolution even in the presence of large amplitude fluctuations of 
    832 the surface height, again so long as $\eta > -H$. 
     777Quite generally, 
     778the time independent range for the \zstar\ coordinate is a very convenient property that 
     779allows for a nearly arbitrary vertical resolution even in 
     780the presence of large amplitude fluctuations of the surface height, again so long as $\eta > -H$. 
    833781%end MOM doc %%% 
    834782 
    835 \newpage 
    836  
    837 % ------------------------------------------------------------------------------------------------------------- 
    838 % Terrain following  coordinate System 
    839 % ------------------------------------------------------------------------------------------------------------- 
     783%% ================================================================================================= 
    840784\subsection{Curvilinear terrain-following \textit{s}--coordinate} 
    841 \label{subsec:PE_sco} 
    842  
    843 % ------------------------------------------------------------------------------------------------------------- 
    844 % Introduction 
    845 % ------------------------------------------------------------------------------------------------------------- 
    846 \subsubsection{Introduction} 
     785\label{subsec:MB_sco} 
    847786 
    848787Several important aspects of the ocean circulation are influenced by bottom topography. 
    849 Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, sills and 
    850 channels that strongly constrain the path of water masses, but more subtle effects exist. 
    851 For example, the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 
     788Of course, the most important is that bottom topography determines deep ocean sub-basins, barriers, 
     789sills and channels that strongly constrain the path of water masses, but more subtle effects exist. 
     790For example, 
     791the topographic $\beta$-effect is usually larger than the planetary one along continental slopes. 
    852792Topographic Rossby waves can be excited and can interact with the mean current. 
    853 In the $z$-coordinate system presented in the previous section (\autoref{sec:PE_zco}), 
     793In the $z$-coordinate system presented in the previous section (\autoref{sec:MB_zco}), 
    854794$z$-surfaces are geopotential surfaces. 
    855795The bottom topography is discretised by steps. 
     
    857797large localized depth gradients associated with large localized vertical velocities. 
    858798The response to such a velocity field often leads to numerical dispersion effects. 
    859 One solution to strongly reduce this error is to use a partial step representation of bottom topography instead of 
    860 a full step one \cite{pacanowski.gnanadesikan_MWR98}. 
     799One solution to strongly reduce this error is to use a partial step representation of 
     800bottom topography instead of a full step one \cite{pacanowski.gnanadesikan_MWR98}. 
    861801Another solution is to introduce a terrain-following coordinate system (hereafter $s$-coordinate). 
    862802 
    863 The $s$-coordinate avoids the discretisation error in the depth field since the layers of 
    864 computation are gradually adjusted with depth to the ocean bottom. 
    865 Relatively small topographic features as well as  gentle, large-scale slopes of the sea floor in the deep ocean, 
    866 which would be ignored in typical $z$-model applications with the largest grid spacing at greatest depths, 
     803The $s$-coordinate avoids the discretisation error in the depth field since 
     804the layers of computation are gradually adjusted with depth to the ocean bottom. 
     805Relatively small topographic features as well as 
     806gentle, large-scale slopes of the sea floor in the deep ocean, 
     807which would be ignored in typical $z$-model applications with 
     808the largest grid spacing at greatest depths, 
    867809can easily be represented (with relatively low vertical resolution). 
    868 A terrain-following model (hereafter $s$-model) also facilitates the modelling of the boundary layer flows over 
    869 a large depth range, which in the framework of the $z$-model would require high vertical resolution over 
     810A terrain-following model (hereafter $s$-model) also facilitates 
     811the modelling of the boundary layer flows over a large depth range, 
     812which in the framework of the $z$-model would require high vertical resolution over 
    870813the whole depth range. 
    871 Moreover, with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as 
     814Moreover, 
     815with a $s$-coordinate it is possible, at least in principle, to have the bottom and the sea surface as 
    872816the only boundaries of the domain (no more lateral boundary condition to specify). 
    873 Nevertheless, a $s$-coordinate also has its drawbacks. Perfectly adapted to a homogeneous ocean, 
     817Nevertheless, a $s$-coordinate also has its drawbacks. 
     818Perfectly adapted to a homogeneous ocean, 
    874819it has strong limitations as soon as stratification is introduced. 
    875820The main two problems come from the truncation error in the horizontal pressure gradient and 
    876821a possibly increased diapycnal diffusion. 
    877 The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx:A}), 
    878  
     822The horizontal pressure force in $s$-coordinate consists of two terms (see \autoref{apdx:SCOORD}), 
    879823\begin{equation} 
    880   \label{eq:PE_p_sco} 
     824  \label{eq:MB_p_sco} 
    881825  \nabla p |_z = \nabla p |_s - \frac{1}{e_3} \pd[p]{s} \nabla z |_s 
    882826\end{equation} 
    883827 
    884 The second term in \autoref{eq:PE_p_sco} depends on the tilt of the coordinate surface and 
     828The second term in \autoref{eq:MB_p_sco} depends on the tilt of the coordinate surface and 
    885829leads to a truncation error that is not present in a $z$-model. 
    886 In the special case of a $\sigma$-coordinate (i.e. a depth-normalised coordinate system $\sigma = z/H$), 
    887 \citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of the magnitude of this truncation error. 
    888 It depends on topographic slope, stratification, horizontal and vertical resolution, the equation of state, 
    889 and the finite difference scheme. 
     830In the special case of a $\sigma$-coordinate 
     831(i.e. a depth-normalised coordinate system $\sigma = z/H$), 
     832\citet{haney_JPO91} and \citet{beckmann.haidvogel_JPO93} have given estimates of 
     833the magnitude of this truncation error. 
     834It depends on topographic slope, stratification, horizontal and vertical resolution, 
     835the equation of state, and the finite difference scheme. 
    890836This error limits the possible topographic slopes that a model can handle at 
    891837a given horizontal and vertical resolution. 
    892838This is a severe restriction for large-scale applications using realistic bottom topography. 
    893 The large-scale slopes require high horizontal resolution, and the computational cost becomes prohibitive. 
     839The large-scale slopes require high horizontal resolution, 
     840and the computational cost becomes prohibitive. 
    894841This problem can be at least partially overcome by mixing $s$-coordinate and 
    895 step-like representation of bottom topography \citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}. 
     842step-like representation of bottom topography 
     843\citep{gerdes_JGR93*a,gerdes_JGR93*b,madec.delecluse.ea_JPO96}. 
    896844However, the definition of the model domain vertical coordinate becomes then a non-trivial thing for 
    897845a realistic bottom topography: 
    898 an envelope topography is defined in $s$-coordinate on which a full or 
    899 partial step bottom topography is then applied in order to adjust the model depth to the observed one 
    900 (see \autoref{sec:DOM_zgr}. 
    901  
    902 For numerical reasons a minimum of diffusion is required along the coordinate surfaces of 
    903 any finite difference model. 
     846an envelope topography is defined in $s$-coordinate on which 
     847a full or partial step bottom topography is then applied in order to 
     848adjust the model depth to the observed one (see \autoref{subsec:DOM_zgr}). 
     849 
     850For numerical reasons a minimum of diffusion is required along 
     851the coordinate surfaces of any finite difference model. 
    904852It causes spurious diapycnal mixing when coordinate surfaces do not coincide with isoneutral surfaces. 
    905853This is the case for a $z$-model as well as for a $s$-model. 
    906 However, density varies more strongly on $s$-surfaces than on horizontal surfaces in regions of 
    907 large topographic slopes, implying larger diapycnal diffusion in a $s$-model than in a $z$-model. 
    908 Whereas such a diapycnal diffusion in a $z$-model tends to weaken horizontal density (pressure) gradients and thus 
    909 the horizontal circulation, it usually reinforces these gradients in a $s$-model, creating spurious circulation. 
    910 For example, imagine an isolated bump of topography in an ocean at rest with a horizontally uniform stratification. 
     854However, density varies more strongly on $s$-surfaces than on horizontal surfaces in 
     855regions of large topographic slopes, 
     856implying larger diapycnal diffusion in a $s$-model than in a $z$-model. 
     857Whereas such a diapycnal diffusion in a $z$-model tends to 
     858weaken horizontal density (pressure) gradients and thus the horizontal circulation, 
     859it usually reinforces these gradients in a $s$-model, creating spurious circulation. 
     860For example, imagine an isolated bump of topography in 
     861an ocean at rest with a horizontally uniform stratification. 
    911862Spurious diffusion along $s$-surfaces will induce a bump of isoneutral surfaces over the topography, 
    912863and thus will generate there a baroclinic eddy. 
    913864In contrast, the ocean will stay at rest in a $z$-model. 
    914 As for the truncation error, the problem can be reduced by introducing the terrain-following coordinate below 
    915 the strongly stratified portion of the water column (\ie\ the main thermocline) \citep{madec.delecluse.ea_JPO96}. 
    916 An alternate solution consists of rotating the lateral diffusive tensor to geopotential or to isoneutral surfaces 
    917 (see \autoref{subsec:PE_ldf}). 
     865As for the truncation error, the problem can be reduced by 
     866introducing the terrain-following coordinate below the strongly stratified portion of the water column 
     867(\ie\ the main thermocline) \citep{madec.delecluse.ea_JPO96}. 
     868An alternate solution consists of rotating the lateral diffusive tensor to 
     869geopotential or to isoneutral surfaces (see \autoref{subsec:MB_ldf}). 
    918870Unfortunately, the slope of isoneutral surfaces relative to the $s$-surfaces can very large, 
    919 strongly exceeding the stability limit of such a operator when it is discretized (see \autoref{chap:LDF}). 
    920  
    921 The $s$-coordinates introduced here \citep{lott.madec.ea_OM90,madec.delecluse.ea_JPO96} differ mainly in two aspects from 
    922 similar models: 
    923 it allows a representation of bottom topography with mixed full or partial step-like/terrain following topography; 
     871strongly exceeding the stability limit of such a operator when it is discretized 
     872(see \autoref{chap:LDF}). 
     873 
     874The $s$-coordinates introduced here \citep{lott.madec.ea_OM90,madec.delecluse.ea_JPO96} 
     875differ mainly in two aspects from similar models: 
     876it allows a representation of bottom topography with 
     877mixed full or partial step-like/terrain following topography; 
    924878It also offers a completely general transformation, $s=s(i,j,z)$ for the vertical coordinate. 
    925879 
    926 % ------------------------------------------------------------------------------------------------------------- 
    927 % Curvilinear z-tilde coordinate System 
    928 % ------------------------------------------------------------------------------------------------------------- 
    929 \subsection{\texorpdfstring{Curvilinear \ztilde-coordinate}{}} 
    930 \label{subsec:PE_zco_tilde} 
    931  
    932 The \ztilde -coordinate has been developed by \citet{leclair.madec_OM11}. 
     880%% ================================================================================================= 
     881\subsection{Curvilinear \ztilde-coordinate} 
     882\label{subsec:MB_zco_tilde} 
     883 
     884The \ztilde-coordinate has been developed by \citet{leclair.madec_OM11}. 
    933885It is available in \NEMO\ since the version 3.4 and is more robust in version 4.0 than previously. 
    934886Nevertheless, it is currently not robust enough to be used in all possible configurations. 
    935887Its use is therefore not recommended. 
    936888 
    937 \newpage 
    938  
    939 % ================================================================ 
    940 % Subgrid Scale Physics 
    941 % ================================================================ 
     889%% ================================================================================================= 
    942890\section{Subgrid scale physics} 
    943 \label{sec:PE_zdf_ldf} 
    944  
    945 The hydrostatic primitive equations describe the behaviour of a geophysical fluid at space and time scales larger than 
    946 a few kilometres in the horizontal, a few meters in the vertical and a few minutes. 
    947 They are usually solved at larger scales: the specified grid spacing and time step of the numerical model. 
     891\label{sec:MB_zdf_ldf} 
     892 
     893The hydrostatic primitive equations describe the behaviour of a geophysical fluid at 
     894space and time scales larger than a few kilometres in the horizontal, 
     895a few meters in the vertical and a few minutes. 
     896They are usually solved at larger scales: the specified grid spacing and time step of 
     897the numerical model. 
    948898The effects of smaller scale motions (coming from the advective terms in the Navier-Stokes equations) 
    949899must be represented entirely in terms of large-scale patterns to close the equations. 
     
    956906small scale processes \textit{in fine} balance the surface input of kinetic energy and heat. 
    957907 
    958 The control exerted by gravity on the flow induces a strong anisotropy between the lateral and vertical motions. 
     908The control exerted by gravity on the flow induces a strong anisotropy between 
     909the lateral and vertical motions. 
    959910Therefore subgrid-scale physics \textbf{D}$^{\vect U}$, $D^{S}$ and $D^{T}$  in 
    960 \autoref{eq:PE_dyn}, \autoref{eq:PE_tra_T} and \autoref{eq:PE_tra_S} are divided into 
    961 a lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 
     911\autoref{eq:MB_PE_dyn}, \autoref{eq:MB_PE_tra_T} and \autoref{eq:MB_PE_tra_S} are divided into 
     912a  lateral part \textbf{D}$^{l \vect U}$, $D^{l S}$ and $D^{l T}$ and 
    962913a vertical part \textbf{D}$^{v \vect U}$, $D^{v S}$ and $D^{v T}$. 
    963 The formulation of these terms and their underlying physics are briefly discussed in the next two subsections. 
    964  
    965 % ------------------------------------------------------------------------------------------------------------- 
    966 % Vertical Subgrid Scale Physics 
    967 % ------------------------------------------------------------------------------------------------------------- 
     914The formulation of these terms and their underlying physics are briefly discussed in 
     915the next two subsections. 
     916 
     917%% ================================================================================================= 
    968918\subsection{Vertical subgrid scale physics} 
    969 \label{subsec:PE_zdf} 
    970  
    971 The model resolution is always larger than the scale at which the major sources of vertical turbulence occur 
    972 (shear instability, internal wave breaking...). 
     919\label{subsec:MB_zdf} 
     920 
     921The model resolution is always larger than the scale at which 
     922the major sources of vertical turbulence occur (shear instability, internal wave breaking...). 
    973923Turbulent motions are thus never explicitly solved, even partially, but always parameterized. 
    974 The vertical turbulent fluxes are assumed to depend linearly on the gradients of large-scale quantities 
    975 (for example, the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$, 
     924The vertical turbulent fluxes are assumed to depend linearly on 
     925the gradients of large-scale quantities (for example, 
     926the turbulent heat flux is given by $\overline{T' w'} = -A^{v T} \partial_z \overline T$, 
    976927where $A^{v T}$ is an eddy coefficient). 
    977928This formulation is analogous to that of molecular diffusion and dissipation. 
     
    981932The resulting vertical momentum and tracer diffusive operators are of second order: 
    982933\begin{equation} 
    983   \label{eq:PE_zdf} 
     934  \label{eq:MB_zdf} 
    984935  \begin{gathered} 
    985     \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \ , \\ 
    986           D^{vT}       = \pd[]{z} \lt( A^{vT} \pd[T]{z}         \rt) \quad \text{and} \quad 
     936    \vect D^{v \vect U} = \pd[]{z} \lt( A^{vm} \pd[\vect U_h]{z} \rt) \text{,} \ 
     937          D^{vT}       = \pd[]{z} \lt( A^{vT} \pd[T]{z}         \rt) \ \text{and} \ 
    987938          D^{vS}       = \pd[]{z} \lt( A^{vT} \pd[S]{z}         \rt) 
    988939  \end{gathered} 
    989940\end{equation} 
    990 where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, respectively. 
     941where $A^{vm}$ and $A^{vT}$ are the vertical eddy viscosity and diffusivity coefficients, 
     942respectively. 
    991943At the sea surface and at the bottom, turbulent fluxes of momentum, heat and salt must be specified 
    992944(see \autoref{chap:SBC} and \autoref{chap:ZDF} and \autoref{sec:TRA_bbl}). 
    993945All the vertical physics is embedded in the specification of the eddy coefficients. 
    994946They can be assumed to be either constant, or function of the local fluid properties 
    995 (\eg\ Richardson number, Brunt-Vais\"{a}l\"{a} frequency, distance from the boundary ...), 
     947(\eg\ Richardson number, Brunt-V\"{a}is\"{a}l\"{a} frequency, distance from the boundary ...), 
    996948or computed from a turbulent closure model. 
    997949The choices available in \NEMO\ are discussed in \autoref{chap:ZDF}). 
    998950 
    999 % ------------------------------------------------------------------------------------------------------------- 
    1000 % Lateral Diffusive and Viscous Operators Formulation 
    1001 % ------------------------------------------------------------------------------------------------------------- 
     951%% ================================================================================================= 
    1002952\subsection{Formulation of the lateral diffusive and viscous operators} 
    1003 \label{subsec:PE_ldf} 
     953\label{subsec:MB_ldf} 
    1004954 
    1005955Lateral turbulence can be roughly divided into a mesoscale turbulence associated with eddies 
    1006956(which can be solved explicitly if the resolution is sufficient since 
    1007957their underlying physics are included in the primitive equations), 
    1008 and a sub mesoscale turbulence which is never explicitly solved even partially, but always parameterized. 
    1009 The formulation of lateral eddy fluxes depends on whether the mesoscale is below or above the grid-spacing 
    1010 (\ie\ the model is eddy-resolving or not). 
     958and a sub mesoscale turbulence which is never explicitly solved even partially, 
     959but always parameterized. 
     960The formulation of lateral eddy fluxes depends on whether 
     961the mesoscale is below or above the grid-spacing (\ie\ the model is eddy-resolving or not). 
    1011962 
    1012963In non-eddy-resolving configurations, the closure is similar to that used for the vertical physics. 
    1013 The lateral turbulent fluxes are assumed to depend linearly on the lateral gradients of large-scale quantities. 
     964The lateral turbulent fluxes are assumed to depend linearly on 
     965the lateral gradients of large-scale quantities. 
    1014966The resulting lateral diffusive and dissipative operators are of second order. 
    1015 Observations show that lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces 
     967Observations show that 
     968lateral mixing induced by mesoscale turbulence tends to be along isopycnal surfaces 
    1016969(or more precisely neutral surfaces \cite{mcdougall_JPO87}) rather than across them. 
    1017 As the slope of neutral surfaces is small in the ocean, a common approximation is to assume that 
    1018 the `lateral' direction is the horizontal, \ie\ the lateral mixing is performed along geopotential surfaces. 
     970As the slope of neutral surfaces is small in the ocean, 
     971a common approximation is to assume that the ``lateral'' direction is the horizontal, 
     972\ie\ the lateral mixing is performed along geopotential surfaces. 
    1019973This leads to a geopotential second order operator for lateral subgrid scale physics. 
    1020 This assumption can be relaxed: the eddy-induced turbulent fluxes can be better approached by assuming that 
     974This assumption can be relaxed: 
     975the eddy-induced turbulent fluxes can be better approached by assuming that 
    1021976they depend linearly on the gradients of large-scale quantities computed along neutral surfaces. 
    1022977In such a case, the diffusive operator is an isoneutral second order operator and 
    1023978it has components in the three space directions. 
    1024 However, 
    1025 both horizontal and isoneutral operators have no effect on mean (\ie\ large scale) potential energy whereas 
     979However, both horizontal and isoneutral operators have no effect on 
     980mean (\ie\ large scale) potential energy whereas 
    1026981potential energy is a main source of turbulence (through baroclinic instabilities). 
    1027982\citet{gent.mcwilliams_JPO90} proposed a parameterisation of mesoscale eddy-induced turbulence which 
    1028983associates an eddy-induced velocity to the isoneutral diffusion. 
    1029984Its mean effect is to reduce the mean potential energy of the ocean. 
    1030 This leads to a formulation of lateral subgrid-scale physics made up of an isoneutral second order operator and 
    1031 an eddy induced advective part. 
     985This leads to a formulation of lateral subgrid-scale physics made up of 
     986an isoneutral second order operator and an eddy induced advective part. 
    1032987In all these lateral diffusive formulations, 
    1033988the specification of the lateral eddy coefficients remains the problematic point as 
    1034 there is no really satisfactory formulation of these coefficients as a function of large-scale features. 
     989there is no really satisfactory formulation of these coefficients as 
     990a function of large-scale features. 
    1035991 
    1036992In eddy-resolving configurations, a second order operator can be used, 
     
    1041997not interfering with the resolved mesoscale activity. 
    1042998Another approach is becoming more and more popular: 
    1043 instead of specifying explicitly a sub-grid scale term in the momentum and tracer time evolution equations, 
     999instead of specifying explicitly a sub-grid scale term in 
     1000the momentum and tracer time evolution equations, 
    10441001one uses an advective scheme which is diffusive enough to maintain the model stability. 
    1045 It must be emphasised that then, all the sub-grid scale physics is included in the formulation of 
    1046 the advection scheme. 
     1002It must be emphasised that then, 
     1003all the sub-grid scale physics is included in the formulation of the advection scheme. 
    10471004 
    10481005All these parameterisations of subgrid scale physics have advantages and drawbacks. 
    1049 They are not all available in \NEMO. For active tracers (temperature and salinity) the main ones are: 
     1006They are not all available in \NEMO. 
     1007For active tracers (temperature and salinity) the main ones are: 
    10501008Laplacian and bilaplacian operators acting along geopotential or iso-neutral surfaces, 
    10511009\citet{gent.mcwilliams_JPO90} parameterisation, and various slightly diffusive advection schemes. 
    1052 For momentum, the main ones are: Laplacian and bilaplacian operators acting along geopotential surfaces, 
     1010For momentum, the main ones are: 
     1011Laplacian and bilaplacian operators acting along geopotential surfaces, 
    10531012and UBS advection schemes when flux form is chosen for the momentum advection. 
    10541013 
     1014%% ================================================================================================= 
    10551015\subsubsection{Lateral laplacian tracer diffusive operator} 
    10561016 
    1057 The lateral Laplacian tracer diffusive operator is defined by (see \autoref{apdx:B}): 
     1017The lateral Laplacian tracer diffusive operator is defined by (see \autoref{apdx:DIFFOPERS}): 
    10581018\begin{equation} 
    1059   \label{eq:PE_iso_tensor} 
    1060   D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad 
    1061   \Re = 
    1062     \begin{pmatrix} 
    1063       1    & 0    & -r_1          \\ 
    1064       0    & 1    & -r_2          \\ 
    1065       -r_1 & -r_2 & r_1^2 + r_2^2 \\ 
    1066     \end{pmatrix} 
     1019  \label{eq:MB_iso_tensor} 
     1020  D^{lT} = \nabla \vect . \lt( A^{lT} \; \Re \; \nabla T \rt) \quad \text{with} \quad \Re = 
     1021  \begin{pmatrix} 
     1022    1    & 0    & -r_1          \\ 
     1023    0    & 1    & -r_2          \\ 
     1024    -r_1 & -r_2 & r_1^2 + r_2^2 \\ 
     1025  \end{pmatrix} 
    10671026\end{equation} 
    10681027where $r_1$ and $r_2$ are the slopes between the surface along which the diffusive operator acts and 
    10691028the model level (\eg\ $z$- or $s$-surfaces). 
    1070 Note that the formulation \autoref{eq:PE_iso_tensor} is exact for 
     1029Note that the formulation \autoref{eq:MB_iso_tensor} is exact for 
    10711030the rotation between geopotential and $s$-surfaces, 
    10721031while it is only an approximation for the rotation between isoneutral and $z$- or $s$-surfaces. 
    1073 Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:PE_iso_tensor} \citep{cox_OM87}. 
    1074 First, the horizontal contribution of the dianeutral mixing is neglected since the ratio between iso and 
    1075 dia-neutral diffusive coefficients is known to be several orders of magnitude smaller than unity. 
     1032Indeed, in the latter case, two assumptions are made to simplify \autoref{eq:MB_iso_tensor} 
     1033\citep{cox_OM87}. 
     1034First, the horizontal contribution of the dianeutral mixing is neglected since 
     1035the ratio between iso and dia-neutral diffusive coefficients is known to be 
     1036several orders of magnitude smaller than unity. 
    10761037Second, the two isoneutral directions of diffusion are assumed to be independent since 
    1077 the slopes are generally less than $10^{-2}$ in the ocean (see \autoref{apdx:B}). 
     1038the slopes are generally less than $10^{-2}$ in the ocean (see \autoref{apdx:DIFFOPERS}). 
    10781039 
    10791040For \textit{iso-level} diffusion, $r_1$ and $r_2 $ are zero. 
     
    10821043For \textit{geopotential} diffusion, 
    10831044$r_1$ and $r_2 $ are the slopes between the geopotential and computational surfaces: 
    1084 they are equal to $\sigma_1$ and $\sigma_2$, respectively (see \autoref{eq:PE_sco_slope}). 
    1085  
    1086 For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between the isoneutral and computational surfaces. 
     1045they are equal to $\sigma_1$ and $\sigma_2$, respectively (see \autoref{eq:MB_sco_slope}). 
     1046 
     1047For \textit{isoneutral} diffusion $r_1$ and $r_2$ are the slopes between 
     1048the isoneutral and computational surfaces. 
    10871049Therefore, they are different quantities, but have similar expressions in $z$- and $s$-coordinates. 
    10881050In $z$-coordinates: 
    10891051\begin{equation} 
    1090   \label{eq:PE_iso_slopes} 
    1091   r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
    1092   r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
     1052  \label{eq:MB_iso_slopes} 
     1053  r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
     1054  r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
    10931055\end{equation} 
    10941056while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 
    10951057 
     1058%% ================================================================================================= 
    10961059\subsubsection{Eddy induced velocity} 
    10971060 
     
    10991062an additional tracer advection is introduced in combination with the isoneutral diffusion of tracers: 
    11001063\[ 
    1101   % \label{eq:PE_iso+eiv} 
     1064  % \label{eq:MB_iso+eiv} 
    11021065  D^{lT} = \nabla \cdot \lt( A^{lT} \; \Re \; \nabla T \rt) + \nabla \cdot \lt( \vect U^\ast \, T \rt) 
    11031066\] 
    11041067where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent, 
    11051068eddy-induced transport velocity. This velocity field is defined by: 
    1106 \begin{gather} 
    1107   % \label{eq:PE_eiv} 
    1108   u^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_1 \rt) \\ 
    1109   v^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \\ 
     1069\[ 
     1070  % \label{eq:MB_eiv} 
     1071  u^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_1 \rt) \quad 
     1072  v^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \quad 
    11101073  w^\ast = - \frac{1}{e_1 e_2} \lt[   \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) 
    11111074                                     + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] 
    1112 \end{gather} 
     1075\] 
    11131076where $A^{eiv}$ is the eddy induced velocity coefficient 
    11141077(or equivalently the isoneutral thickness diffusivity coefficient), 
    1115 and $\tilde r_1$ and $\tilde r_2$ are the slopes between isoneutral and \textit{geopotential} surfaces. 
    1116 Their values are thus independent of the vertical coordinate, but their expression depends on the coordinate: 
    1117 \begin{align} 
    1118   \label{eq:PE_slopes_eiv} 
     1078and $\tilde r_1$ and $\tilde r_2$ are the slopes between 
     1079isoneutral and \textit{geopotential} surfaces. 
     1080Their values are thus independent of the vertical coordinate, 
     1081but their expression depends on the coordinate: 
     1082\begin{equation} 
     1083  \label{eq:MB_slopes_eiv} 
    11191084  \tilde{r}_n = 
    11201085    \begin{cases} 
     
    11231088    \end{cases} 
    11241089  \quad \text{where~} n = 1, 2 
    1125 \end{align} 
     1090\end{equation} 
    11261091 
    11271092The normal component of the eddy induced velocity is zero at all the boundaries. 
    1128 This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in the vicinity of 
    1129 the boundaries. 
     1093This can be achieved in a model by tapering either the eddy coefficient or the slopes to zero in 
     1094the vicinity of the boundaries. 
    11301095The latter strategy is used in \NEMO\ (cf. \autoref{chap:LDF}). 
    11311096 
     1097%% ================================================================================================= 
    11321098\subsubsection{Lateral bilaplacian tracer diffusive operator} 
    11331099 
    11341100The lateral bilaplacian tracer diffusive operator is defined by: 
    11351101\[ 
    1136   % \label{eq:PE_bilapT} 
     1102  % \label{eq:MB_bilapT} 
    11371103  D^{lT}= - \Delta \; (\Delta T) \quad \text{where} \quad 
    11381104  \Delta \bullet = \nabla \lt( \sqrt{B^{lT}} \; \Re \; \nabla \bullet \rt) 
    11391105\] 
    1140 It is the Laplacian operator given by \autoref{eq:PE_iso_tensor} applied twice with 
     1106It is the Laplacian operator given by \autoref{eq:MB_iso_tensor} applied twice with 
    11411107the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 
    11421108 
     1109%% ================================================================================================= 
    11431110\subsubsection{Lateral Laplacian momentum diffusive operator} 
    11441111 
    11451112The Laplacian momentum diffusive operator along $z$- or $s$-surfaces is found by 
    1146 applying \autoref{eq:PE_lap_vector} to the horizontal velocity vector (see \autoref{apdx:B}): 
     1113applying \autoref{eq:MB_lap_vector} to the horizontal velocity vector (see \autoref{apdx:DIFFOPERS}): 
    11471114\begin{align*} 
    1148   % \label{eq:PE_lapU} 
     1115  % \label{eq:MB_lapU} 
    11491116  \vect D^{l \vect U} &=   \nabla_h        \big( A^{lm}    \chi             \big) 
    11501117                         - \nabla_h \times \big( A^{lm} \, \zeta \; \vect k \big) \\ 
     
    11551122\end{align*} 
    11561123 
    1157 Such a formulation ensures a complete separation between the vorticity and horizontal divergence fields 
    1158 (see \autoref{apdx:C}). 
     1124Such a formulation ensures a complete separation between 
     1125the vorticity and horizontal divergence fields (see \autoref{apdx:INVARIANTS}). 
    11591126Unfortunately, it is only available in \textit{iso-level} direction. 
    11601127When a rotation is required 
    1161 (\ie\ geopotential diffusion in $s$-coordinates or isoneutral diffusion in both $z$- and $s$-coordinates), 
    1162 the $u$ and $v$-fields are considered as independent scalar fields, so that the diffusive operator is given by: 
    1163 \begin{gather*} 
    1164   % \label{eq:PE_lapU_iso} 
    1165     D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\ 
     1128(\ie\ geopotential diffusion in $s$-coordinates or 
     1129isoneutral diffusion in both $z$- and $s$-coordinates), 
     1130the $u$ and $v$-fields are considered as independent scalar fields, 
     1131so that the diffusive operator is given by: 
     1132\[ 
     1133  % \label{eq:MB_lapU_iso} 
     1134    D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \quad 
    11661135    D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) 
    1167 \end{gather*} 
    1168 where $\Re$ is given by \autoref{eq:PE_iso_tensor}. 
     1136\] 
     1137where $\Re$ is given by \autoref{eq:MB_iso_tensor}. 
    11691138It is the same expression as those used for diffusive operator on tracers. 
    11701139It must be emphasised that such a formulation is only exact in a Cartesian coordinate system, 
     
    11731142a geographical coordinate system \citep{lengaigne.madec.ea_JGR03}. 
    11741143 
     1144%% ================================================================================================= 
    11751145\subsubsection{Lateral bilaplacian momentum diffusive operator} 
    11761146 
    1177 As for tracers, the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with 
     1147As for tracers, 
     1148the bilaplacian order momentum diffusive operator is a re-entering Laplacian operator with 
    11781149the harmonic eddy diffusion coefficient set to the square root of the biharmonic one. 
    11791150Nevertheless it is currently not available in the iso-neutral case. 
    11801151 
    1181 \biblio 
    1182  
    1183 \pindex 
     1152\onlyinsubfile{\input{../../global/epilogue}} 
    11841153 
    11851154\end{document} 
Note: See TracChangeset for help on using the changeset viewer.