Skip to content
Snippets Groups Projects
apdx_DOMAINcfg.tex 42.8 KiB
Newer Older
\documentclass[../main/NEMO_manual]{subfiles}

\begin{document}

\chapter{A brief guide to the DOMAINcfg tool}
\label{apdx:DOMCFG}

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{\textwidth}{l||X|X}
    Release     & Author(s)            & Modifications                                                 \\
    \hline
     {\em  5.0} & {\em Katherine Hutchinson and Pierre Mathiot} & {\em update to NEMO 5.0 } \\
    {\em  4.2} & {\em Pierre Mathiot} & {\em Add ice shelf and closed sea option description        } \\
    {\em   4.0} & {\em  Andrew Coward} & {\em Creation from materials removed from \autoref{chap:DOM}
                                              that are still relevant to the DOMAINcfg tool
                                              when setting up new domains                            }
  \end{tabularx}
}

\clearpage

This appendix briefly describes some of the options available in the
\forcode{DOMAINcfg} tool mentioned in \autoref{chap:DOM}.

Please note that there are plans to update and clean this tool and provide better step-by-step instructions. 
In the interim, for practical guidelines on the use of the tool, please see the \href{https://forge.nemo-ocean.eu/nemo/nemo/-/blob/4.2.2/tools/DOMAINcfg/README.rst}{DOMAINcfg README}. 

The \texttt{DOMAINcfg} tool allows the user to define some
horizontal and vertical grids through additional namelist parameters. Explanations of
these parameters are retained here for reference pending better documentation for
\forcode{DOMAINcfg}. Please note that the namelist blocks named in this appendix refer to
those read by \forcode{DOMAINcfg} via its own \forcode{namelist_ref} and
\forcode{namelist_cfg} files. Although, due to their origins, these namelists share names
with those used by \NEMO, they are not interchangeable and should be considered independent
of those described elsewhere in this manual.

%% =================================================================================================
\section{Choice of horizontal grid}
\label{sec:DOMCFG_hor}

\begin{listing}
  \begin{forlines}
!-----------------------------------------------------------------------
&namdom        !   space and time domain (bathymetry, mesh, timestep)
!-----------------------------------------------------------------------
   ln_read_cfg = .false.   !  Read from a domain_cfg file
   nn_bathy    =    1      ! = 0 compute analyticaly
                           ! = 1 read the bathymetry file
                           ! = 2 compute from external bathymetry
                           ! = 3 compute from parent (if "key_agrif")
   nn_interp   =    1                          ! type of interpolation (nn_bathy =2)
   cn_domcfg   = ' '       ! Name of the domain_cfg input file
   cn_fcoord   =  'coordinates.nc'             ! external coordinates file (jphgr_msh = 0)
   cn_topo     =  'bathy_meter.nc           '  ! external topo file (nn_bathy =1/2)
   cn_topolvl  =  'bathy_level.nc           '  ! external topo file (nn_bathy =1)
   cn_fisfd    =  'isf_draft_meter.nc'         ! external isf draft (nn_bathy =1 and ln_isfcav = .true.)
   cn_bath     =  'Bathymetry'                 ! topo name in file  (nn_bathy =1/2)
   cn_bathlvl  =  'Bathy_level'                ! lvl name in file   (nn_bathy =1)
   cn_visfd    =  'isf_draft'                  ! isf draft variable (nn_bathy =1 and ln_isfcav = .true.)
   cn_lon      =  'nav_lon'                    ! lon  name in file  (nn_bathy =2)
   cn_lat      =  'nav_lat'                    ! lat  name in file  (nn_bathy =2)
   rn_scale    =    1.     !  multiplicative factor to account for possibly negative input bathymetry (agrif only) 
   rn_bathy    =    0.     !  value of the bathymetry. if (=0) bottom flat at jpkm1
   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
   rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
   rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
   rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
                           !
   rn_rdt      = 5760.     !  time step for the dynamics (and tracer if nn_acc=0)
   rn_atfp     =    0.1    !  asselin time filter parameter
   ln_crs      = .false.      !  Logical switch for coarsening module
   jphgr_msh   =       0               !  type of horizontal mesh
                                       !  = 0 curvilinear coordinate on the sphere read in coordinate.nc
                                       !  = 1 geographical mesh on the sphere with regular grid-spacing
                                       !  = 2 f-plane with regular grid-spacing
                                       !  = 3 beta-plane with regular grid-spacing
                                       !  = 4 Mercator grid with T/U point at the equator
   ppglam0     =       0.0             !  longitude of first raw and column T-point (jphgr_msh = 1)
   ppgphi0     =     -35.0             ! latitude  of first raw and column T-point (jphgr_msh = 1)
   ppe1_deg    =       1.0             !  zonal      grid-spacing (degrees)
   ppe2_deg    =       0.5             !  meridional grid-spacing (degrees)
   ppe1_m      =    5000.0             !  zonal      grid-spacing (degrees)
   ppe2_m      =    5000.0             !  meridional grid-spacing (degrees)
   ppsur       =    -4762.96143546300  !  ORCA r4, r2 and r05 coefficients
   ppa0        =      255.58049070440  ! (default coefficients)
   ppa1        =      245.58132232490  !
   ppkth       =       21.43336197938  !
   ppacr       =        3.0            !
   ppdzmin     =       10.             !  Minimum vertical spacing
   pphmax      =     5000.             !  Maximum depth
   ldbletanh   =    .TRUE.             !  Use/do not use double tanf function for vertical coordinates
   ppa2        =      100.760928500000 !  Double tanh function parameters
   ppkth2      =       48.029893720000 !
   ppacr2      =       13.000000000000 !
/
  \end{forlines}
  \caption{\forcode{&namdom_domcfg}}
  \label{lst:namdom_domcfg}
\end{listing}

The user has three options available in defining a horizontal grid, which involve the
namelist variable \np{jphgr_msh}{jphgr\_msh} (\autoref{lst:namdom_domcfg}).

 \item [{\np{jphgr_msh}{jphgr\_msh}=0}]  The most general curvilinear orthogonal grids.
  The coordinates and their first derivatives with respect to $i$ and $j$ are provided
  in a input file (\textit{coordinates.nc}), read in \rou{hgr\_read} subroutine of the domhgr module.
  This is now the only option available within \NEMO\ itself from v4.0 onwards.
\item [{\np{jphgr_msh}{jphgr\_msh}=1 to 4}] A few simple analytical grids are provided (see below).
  For other analytical grids, the \mdl{domhgr} module (\texttt{DOMAINcfg} variant) must be
  modified by the user. In most cases, modifying the \mdl{usrdef\_hgr} module of \NEMO\ is
  a better alternative since this is designed to allow simple analytical domains to be
  configured and used without the need for external data files.
\end{description}

There are two simple cases of geographical grids on the sphere. With
\np[=1]{jphgr_msh}{jphgr\_msh}, the grid (expressed in degrees) is regular in space,
with grid sizes specified by parameters \np{ppe1_deg}{ppe1\_deg} and \np{ppe2_deg}{ppe2\_deg},
respectively. Such a geographical grid can be very anisotropic at high latitudes
because of the convergence of meridians (the zonal scale factors $e_1$
become much smaller than the meridional scale factors $e_2$). The Mercator
grid (\np{jphgr_msh}{jphgr\_msh}=4) avoids this anisotropy by refining the meridional scale
factors in the same way as the zonal ones. In this case, meridional scale factors
and latitudes are calculated analytically using the formulae appropriate for
a Mercator projection, based on \np{ppe1_deg}{ppe1\_deg} which is a reference grid spacing
at the equator (this applies even when the geographical equator is situated outside
the model domain).

In these two cases (\np{jphgr_msh}{jphgr\_msh}=1 or 4), the grid position is defined by the
longitude and latitude of the south-westernmost point (\np{ppglam0}{ppglam0}
and \np{ppgphi0}{ppgphi0}). Note that for the Mercator grid the user needs only to provide
an approximate starting latitude: the real latitude will be recalculated analytically,
in order to ensure that the equator corresponds to line passing through $t$-
and $u$-points.

Rectangular grids ignoring the spherical geometry are defined with
\np{jphgr_msh}{jphgr\_msh} = 2 and 3. The domain is either an $f$-plane (\np{jphgr_msh}{jphgr\_msh} = 2,
Coriolis factor is constant) or a beta-plane (\np{jphgr_msh}{jphgr\_msh} = 3, the Coriolis factor
is linear in the $j$-direction). The grid size is uniform in meters in each direction,
and given by the parameters \np{ppe1_m}{ppe1\_m} and \np{ppe2_m}{ppe2\_m} respectively.
The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero
with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers,
and the second $t$-point corresponds to coordinate $gphit=0$. The input
variable \np{ppglam0}{ppglam0} is ignored. \np{ppgphi0}{ppgphi0} is used to set the reference
latitude for computation of the Coriolis parameter. In the case of the beta plane,
\np{ppgphi0}{ppgphi0} corresponds to the center of the domain. 

%% =================================================================================================
\section{Vertical grid}
\label{sec:DOMCFG_vert}

%% =================================================================================================
\subsection{Vertical reference coordinate}
\label{sec:DOMCFG_zref}

\begin{figure}[!tb]
  \centering
  \includegraphics[width=0.66\textwidth]{DOMCFG_zgr}
  \caption[DOMAINcfg: default vertical mesh for ORCA2]{
    Default vertical mesh for ORCA2: 30 ocean levels (L30).
    Vertical level functions for (a) T-point depth and (b) the associated scale factor for
    the $z$-coordinate case.}
  \label{fig:DOMCFG_zgr}
\end{figure}

The reference coordinate transformation $z_0(k)$ defines the arrays $gdept\_0$ and
$gdepw\_0$ for $t$- and $w$-points, respectively. See \autoref{sec:DOMCFG_sco} for the
S-coordinate options.  As indicated on \autoref{fig:DOM_index_vert} \texttt{jpk} is the number of
$w$-levels.  $gdepw\_0(1)$ is the ocean surface.  There are at most \texttt{jpk}-1 $t$-points
inside the ocean, the additional $t$-point at $jk = jpk$ is below the sea floor and is not

Since version 4.0, the default definition of cell depths for both $w$- and $t$-levels is done in a discrete sense
where $e_3^0T = d_k(z^W_0)$ and $e_3^0W = d_k(z^T_0)$, by setting \np[=.true.]{ln_e3_dep}{ln\_e3\_dep}
in the (\autoref{lst:namdom_domcfg}) (\texttt{DOMAINcfg} \texttt{\&namdom} variant).
(namely $z^T_0(k) = 0.5 (z^W_0(k)+z^W_0(k+1))$ and $ e_3^0W = d_k(z^T_0)$).
For backward compatibility with version 3.6, by setting \np[=.false.]{ln_e3_dep}{ln\_e3\_dep}
the vertical location of $w$- and $t$-levels are defined using the former analytic
expression of the depth $z_0(k)$, whose analytical derivative with respect to $k$ provides
the vertical scale factors.
All the above operations are done through statement functions in the routine \mdl{domzgr} (\texttt{DOMAINcfg} variant).

It is possible to define a simple regular vertical grid by giving zero stretching
(\np[=0]{ppacr}{ppacr}).  In that case, the parameters \texttt{jpk} (number of $w$-levels)
and \np{pphmax}{pphmax} (maximum ocean depth in meters) fully define the grid.

For climate-related studies it is often desirable to concentrate the vertical resolution
near the ocean surface.  The following function is proposed as a standard for a
$z$-coordinate (with either full or partial steps):
\begin{gather}
  \label{eq:DOMCFG_zgr_ana_1}
    z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; h_{cr} \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\
    e_3^0(k) = \lt|    - h_0      -    h_1 \;           \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt|
where $k = 1$ to \texttt{jpk} for $w$-levels and $k = 1$ to \texttt{jpk-1} for $t-$levels.
Such an expression allows us to define a nearly uniform vertical location of levels at the ocean
top and bottom with a smooth hyperbolic tangent transition in between (\autoref{fig:DOMCFG_zgr}).

A double hyperbolic tangent version (\np[=.true.]{ldbletanh}{ldbletanh}) is also available
which permits finer control and is used, typically, to obtain a well resolved upper ocean
without compromising on resolution at depth using a moderate number of levels.

\begin{gather}
  \label{eq:DOMCFG_zgr_ana_1b}
    \begin{split}
    z_0  (k) = h_{sur} - h_0 \; k &- \; h_1  \; h_{cr} \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\
                             \;   &- \; h2_1 \; h_{cr} \; \log  \big[ \cosh ((k - h2_{th}) / h2_{cr}) \big]
    \end{split}
\end{gather}
\begin{gather}
    \begin{split}
    e_3^0(k) = \big|    - h_0    &-   h_1 \; \tanh \big[       (k - h_{th})  / h_{cr}   \big]  \\
                                 &-  h2_1 \; \tanh \big[       (k - h2_{th}) / h2_{cr}  \big] \big|
    \end{split}
\end{gather}

If the ice shelf cavities are opened (\np[=.true.]{ln_isfcav}{ln\_isfcav}), the definition
of $z_0$ is the same.  However, definition of $e_3^0$ at $t$- and $w$-points is
respectively changed to:
\begin{equation}
  \label{eq:DOMCFG_zgr_ana_2}
  \begin{split}
    e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\
    e_3^W(k) &= z_T (k    ) - z_T (k - 1)
  \end{split}
\end{equation}

This formulation decreases the self-generated circulation into the ice shelf cavity
(which can, in extreme case, leads to numerical instability). This is now the recommended formulation for all configurations using v4.0 onwards. The analytical derivation of thicknesses is maintained for backwards compatibility.

The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface
(bottom) layers and a depth which varies from 0 at the sea surface to a minimum of
$-5000~m$.  This leads to the following conditions:

\begin{equation}
  \label{eq:DOMCFG_zgr_coef}
  \begin{array}{ll}
    e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\
    e_3 (jpk - 1/2) = 500. & z(jpk) = -5000.
  \end{array}
\end{equation}

With the choice of the stretching $h_{cr} = 3$ and the number of levels \texttt{jpk}~$= 31$,
the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in
\autoref{eq:DOMCFG_zgr_ana_2} have been determined such that \autoref{eq:DOMCFG_zgr_coef}
is satisfied, through an optimisation procedure using a bisection method.
For the first standard ORCA2 vertical grid this led to the following values:
$h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$.
The resulting depths and scale factors as a function of the model levels are shown in
\autoref{fig:DOMCFG_zgr} and given in \autoref{tab:DOMCFG_orca_zgr}.
Those values correspond to the parameters \np{ppsur}{ppsur}, \np{ppa0}{ppa0}, \np{ppa1}{ppa1}, \np{ppkth}{ppkth} in the \texttt{DOMAINcfg} variant of \texttt{\&namdom} (\autoref{lst:namdom_domcfg}).

Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to
recalculate them.  In that case the user sets \np{ppsur}{ppsur}~$=$~\np{ppa0}{ppa0}~$=$~\np{ppa1}{ppa1}~$=
999999$., in \nam{cfg}{cfg} namelist, and specifies instead the four following parameters:
\begin{itemize}
\item \np{ppacr}{ppacr}~$= h_{cr}$: stretching factor (nondimensional).
  The larger \np{ppacr}{ppacr}, the smaller the stretching.
  Values from $3$ to $10$ are usual.
\item \np{ppkth}{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs
  (nondimensional, usually of order 1/2 or 2/3 of \texttt{jpk})
\item \np{ppdzmin}{ppdzmin}: minimum thickness for the top layer (in meters).
\item \np{pphmax}{pphmax}: total depth of the ocean (meters).
\end{itemize}

As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are:
\texttt{jpk}~$= 46$, \np{ppacr}{ppacr}~$= 9$, \np{ppkth}{ppkth}~$= 23.563$, \np{ppdzmin}{ppdzmin}~$= 6~m$,
\np{pphmax}{pphmax}~$= 5750~m$.

\begin{table}
  \centering
  \begin{tabular}{c||r|r|r|r}
    \hline
    \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\
    \hline
    1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\
    \hline
    2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\
    \hline
    3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\
    \hline
    4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\
    \hline
    5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\
    \hline
    6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\
    \hline
    7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\
    \hline
    8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\
    \hline
    9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\
    \hline
    10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\
    \hline
    11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\
    \hline
    12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\
    \hline
    13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\
    \hline
    14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\
    \hline
    15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\
    \hline
    16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\
    \hline
    17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\
    \hline
    18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\
    \hline
    19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\
    \hline
    20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\
    \hline
    21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\
    \hline
    22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\
    \hline
    23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\
    \hline
    24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\
    \hline
    25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\
    \hline
    26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\
    \hline
    27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\
    \hline
    28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\
    \hline
    29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\
    \hline
    30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\
    \hline
    31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\
    \hline
  \end{tabular}
  \caption[Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration]{
    Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as
    computed from \autoref{eq:DOMCFG_zgr_ana_2} using
    the coefficients given in \autoref{eq:DOMCFG_zgr_coef}}
  \label{tab:DOMCFG_orca_zgr}
\end{table}
%%%YY
%% % -------------------------------------------------------------------------------------------------------------
%% %        Meter Bathymetry
%% % -------------------------------------------------------------------------------------------------------------
%% =================================================================================================
\subsection{Model bathymetry}
\label{subsec:DOMCFG_bathy}

The bathymetry can either be defined analytically by hand in \texttt{DOMAIncfg} variant of \mdl{domzgr} in the \rou{zgr\_bat} subroutine (\np[ <= 0]{nn_bathy}{nn\_bathy}) with some predefined bathymetry or from a netcdf file provided by the user (\np[ > 0]{nn_bathy}{nn\_bathy}).

For the case \np[ <= 0]{nn_bathy}{nn\_bathy}:
\begin{description}
\item [{\np[=0]{nn_bathy}{nn\_bathy}}]: a flat-bottom domain is defined.
  The total depth $z_w (jpk)$ is given by the coordinate transformation.
  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}{jperio}.
\item [{\np[=-1]{nn_bathy}{nn\_bathy}}]: There are two predefined bathymetry. One for the OVERFLOW, and one for the DOME test cases (see user guide). If the configuration name (\np{cp_cfg}{cp\_cfg} is different from DOME or OVERFLOW, the default is a flat bottom with a random noise. The user is strongly encourage to check \rou{zgr\_bat} to see if it fits exactly its need (many parameter hard coded for each cases).
\end{description}

For the case \np[ > 0]{nn_bathy}{nn\_bathy}:
\begin{description}
\item [{\np[=1]{nn_bathy}{nn\_bathy}}]: There is two cases. For z or $\sigma$ coordinates (\np{ln_zps}{ln\_zps} or \np{ln_sco}{ln\_sco}), the user needs to provide a bathymetry and ice shelf draft depth variables (\np{cn_bath}{cn\_bath} and \np{cn_visfd}{cn\_visfd}) in netcdf files (\np{cn_topo}{cn\_topo} and \np{cn_fisfd}{cn\_fisfd}). Ice shelf draft variable is optional as it depends if \np{ln_isfcav}{ln\_isfcav} is activated or not. Each files need to provide data on each grid point of the 2D model grid. The depth convention for both variables is positive, in meters and the continent are defined by a 0 or negative value.\\
  The bathymetry is usually built by interpolating a standard bathymetry product (\eg\ ETOPO or GEBCO) onto
  the horizontal ocean mesh. For the case full cell (\np{ln_zco}{ln\_zco}), the user need to provide a bottom level variables (\np{cn_bathlvl}{cn\_bathlvl} in a netcdf file (\np{cn_topolvl}{cn\_topolvl}). Representation of the ice shel cavities is not implemented in this case.\\
  In both cases, there is specific hand edits for the ORCA2 configuration (\np[='orca']{cp_cfg}{cp\_cfg} and \np[=2]{jp_cfg}{jp\_cfg}) for Gibraltar and Bab el Mandeb straits.
 \item [{\np[>1]{nn_bathy}{nn\_bathy}}]: ???
\end{description}

%% =================================================================================================
\subsection{Choice of vertical grid}
\label{sec:DOMCFG_vgrd}

After reading the bathymetry, the algorithm for vertical grid definition differs between the different options:
\begin{description}
\item [\forcode{ln_zco = .true.}] set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$ where $z_0(k)$ is the closest match to the depth at $(i,j)$.
\item [\forcode{ln_zps = .true.}] set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level (and shallowest level if \np{ln_isfcav}{ln\_isfcav}) at
  each $(i,j)$ point using the bathymetry (and ice shelf draftif \np{ln_isfcav}{ln\_isfcav}), to obtain the final three-dimensional depth and scale factor arrays.
\item [\forcode{ln_sco = .true.}] smooth the bathymetry to fulfill the hydrostatic consistency criteria and
  set the three-dimensional transformation.
\item [\forcode{s-z and s-zps}] smooth the bathymetry to fulfill the hydrostatic consistency criteria and
  set the three-dimensional transformation $z(i,j,k)$,
  and possibly introduce masking of extra land points to better fit the original bathymetry file.
\end{description}

%% =================================================================================================
\subsubsection[$Z$-coordinate with uniform thickness levels (\forcode{ln_zco})]{$Z$-coordinate with uniform thickness levels (\protect\np{ln_zco}{ln\_zco})}
\label{subsec:DOMCFG_zco}

With this option the model topography can be fully described by the reference vertical
coordinate and a 2D integer field giving the number of wet levels at each location
(\forcode{bathy_level}). The resulting match to the real topography is likely to be poor
though (especially with thick, deep levels) and slopes poorly represented. This option is
rarely used in modern simulations but it can be useful for testing purposes.

%% =================================================================================================
\subsubsection[$Z$-coordinate with partial step (\forcode{ln_zps})]{$Z$-coordinate with partial step (\protect\np{ln_zps}{ln\_zps})}
\label{subsec:DOMCFG_zps}

In $z$-coordinate partial step, the depths of the model levels are defined by the
reference analytical function $z_0(k)$ as described in \autoref{sec:DOMCFG_zref},
\textit{except} in the bottom layer.  The thickness of the bottom layer is allowed to vary
as a function of geographical location $(\lambda,\varphi)$ to allow a better
representation of the bathymetry, especially in the case of small slopes (where the
bathymetry varies by less than one level thickness from one grid point to the next).  The
reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry.
With partial steps, layers from 1 to \texttt{jpk-2} can have a thickness smaller than
$e_{3t}(jk)$.

The model deepest layer (\texttt{jpk-1}) is allowed to have either a smaller or larger
thickness than $e_{3t}(jpk)$: the maximum thickness allowed is $2*e_{3t}(jpk - 1)$.

In case of ice shelf, partial step are allowed at the top interface using the same methodology. 

This has to be kept in mind when specifying values in  \texttt{\&namdom} (\texttt{DOMAIncfg} variant ;  \autoref{lst:namdom_domcfg}), such as the maximum depth \np{pphmax}{pphmax} in partial steps.

For example, with \np{pphmax}{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean
depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being
$250~m$).  Two variables in the namdom namelist are used to define the partial step
vertical grid.  The mimimum water thickness (in meters) allowed for a cell partially
filled with bathymetry at level jk is the minimum of \np{rn_e3zps_min}{rn\_e3zps\_min} (thickness in
meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn_e3zps_rat}{rn\_e3zps\_rat} (a fraction, usually 10\%, of
the default thickness $e_{3t}(jk)$).

%% =================================================================================================
\subsubsection[$S$-coordinate (\forcode{ln_sco})]{$S$-coordinate (\protect\np{ln_sco}{ln\_sco})}
\label{sec:DOMCFG_sco}

\begin{listing}
  \caption{\forcode{&namzgr_sco_domcfg}}
  \label{lst:namzgr_sco_domcfg}
  \begin{forlines}
!-----------------------------------------------------------------------
&namzgr_sco    !   s-coordinate or hybrid z-s-coordinate                (default: OFF)
!-----------------------------------------------------------------------
   ln_s_sh94   = .false.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
                           !  stretching coefficients for all functions
   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
                        !!!!!!!  Envelop bathymetry
   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
   rn_bb       =    0.8    !  stretching with SH94 s-sigma
                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
   rn_zs       =    1.0    !  depth of surface grid box
                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
   rn_zb_b     =   -0.2    !  offset for calculating Zb
                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)
/
  \end{forlines}
\end{listing}

Options are defined in \forcode{&namzgr_sco} (\texttt{DOMAINcfg} only).
In $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}), the depth and thickness of the model levels are defined from
the product of a depth field and either a stretching function or its derivative, respectively:

\begin{align*}
  % \label{eq:DOMCFG_sco_ana}
  z(k)   &= h(i,j) \; z_0 (k) \\
  e_3(k) &= h(i,j) \; z_0'(k)
\end{align*}

where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and
$z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom.
The depth field $h$ is not necessary the ocean depth,
since a mixed step-like and bottom-following representation of the topography can be used
(\autoref{fig:DOM_z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:DOM_z_zps_s_sps}).
The namelist parameter \np{rn_rmax}{rn\_rmax} determines the slope at which
the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.
The coordinate can also be hybridised by specifying \np{rn_sbot_min}{rn\_sbot\_min} and \np{rn_sbot_max}{rn\_sbot\_max} as
the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.

Options for stretching the coordinate are provided as examples,
but care must be taken to ensure that the vertical stretch used is appropriate for the application.

The original default \NEMO\ s-coordinate stretching is available if neither of the other options are specified as true
(\np[=.false.]{ln_s_sh94}{ln\_s\_sh94} and \np[=.false.]{ln_s_SF12}{ln\_s\_SF12}).
This uses a depth independent $\tanh$ function for the stretching \citep{madec.delecluse.ea_JPO96}:

\[
  z = s_{min} + C (s) (H - s_{min})
  % \label{eq:DOMCFG_SH94_1}
\]

where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and
allows a $z$-coordinate to be placed on top of the stretched coordinate,
and $z$ is the depth (negative down from the asea surface).
\begin{gather*}
  s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1
  % \label{eq:DOMCFG_s}
 \\
 \label{eq:DOMCFG_sco_function}
  C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)}
\end{gather*}

A stretching function,
modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np[=.true.]{ln_s_sh94}{ln\_s\_sh94}),
is also available and is more commonly used for shelf seas modelling:

\[
  C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)}
         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)}
                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)}
 \label{eq:DOMCFG_SH94_2}
\]

\begin{figure}[!ht]
  \centering
  \includegraphics[width=0.80\textwidth]{DOMCFG_sco_function}
  \caption[DOMAINcfg: examples of the stretching function applied to a seamount]{
    Examples of the stretching function applied to a seamount;
    from left to right: surface, surface and bottom, and bottom intensified resolutions}
  \label{fig:DOMCFG_sco_function}
\end{figure}

where $H_c$ is the critical depth (\np{rn_hc}{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to
the stretched coordinate, and $\theta$ (\np{rn_theta}{rn\_theta}) and $b$ (\np{rn_bb}{rn\_bb}) are the surface and
bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$.
$b$ has been designed to allow surface and/or bottom increase of the vertical resolution
(\autoref{fig:DOMCFG_sco_function}).

Another example has been provided at version 3.5 (\np{ln_s_sf12}{ln\_s\_sf12}) that allows a fixed surface resolution in
an analytical terrain-following stretching \citet{siddorn.furner_OM13}.
In this case the a stretching function $\gamma$ is defined such that:

\begin{equation}
  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1
  % \label{eq:DOMCFG_z}
\end{equation}

The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:

\begin{gather*}
  % \label{eq:DOMCFG_gamma_deriv}
  \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt)
           + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\
  \intertext{Where:}
 \label{eq:DOMCFG_gamma}
  f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2}
  \quad \text{and} \quad \sigma = \frac{k}{n - 1}
\end{gather*}

This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
the user prescribed stretching parameter $\alpha$ (\np{rn_alpha}{rn\_alpha}) that stretches towards
the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and
user prescribed surface (\np{rn_zs}{rn\_zs}) and bottom depths.
The bottom cell depth in this example is given as a function of water depth:

\[
  % \label{eq:DOMCFG_zb}
  Z_b = h a + b
\]

where the namelist parameters \np{rn_zb_a}{rn\_zb\_a} and \np{rn_zb_b}{rn\_zb\_b} are $a$ and $b$ respectively.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=0.66\textwidth]{DOMCFG_compare_coordinates_surface}
  \caption[DOMAINcfg: comparison of $s$- and $z$-coordinate]{
    A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines),
    a 50 level $Z$-coordinate (contoured surfaces) and
    the \citet{siddorn.furner_OM13} $S$-coordinate (dashed lines) in the surface $100~m$ for
    a idealised bathymetry that goes from $50~m$ to $5500~m$ depth.
    For clarity every third coordinate surface is shown.}
  \label{fig:DOMCFG_fig_compare_coordinates_surface}
\end{figure}
 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>

This gives a smooth analytical stretching in computational space that is constrained to
given specified surface and bottom grid cell thicknesses in real space.
This is not to be confused with the hybrid schemes that
superimpose geopotential coordinates on terrain following coordinates thus
creating a non-analytical vertical coordinate that
therefore may suffer from large gradients in the vertical resolutions.
This stretching is less straightforward to implement than the \citet{song.haidvogel_JCP94} stretching,
but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.

As with the \citet{song.haidvogel_JCP94} stretching the stretch is only applied at depths greater than
the critical depth $h_c$.
In this example two options are available in depths shallower than $h_c$,
with pure sigma being applied if the \np{ln_sigcrit}{ln\_sigcrit} is true and pure z-coordinates if it is false
(the z-coordinate being equal to the depths of the stretched coordinate at $h_c$).

Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as
large slopes lead to hydrostatic consistency.
A hydrostatic consistency parameter diagnostic following \citet{haney_JPO91} has been implemented,
and is output as part of the model mesh file at the start of the run.

%% =================================================================================================
\subsubsection[\zstar- or \sstar-coordinate (\forcode{ln_linssh})]{\zstar- or \sstar-coordinate (\protect\np{ln_linssh}{ln\_linssh})}
\label{subsec:DOMCFG_zgr_star}

This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO\ web site.

\section{Ice shelf cavity definition}
\label{subsec:zgrisf}

\begin{listing}
  \caption{\forcode{&namzgr_isf}}
  \label{lst:namzgr_isf}
  \begin{forlines}
!-----------------------------------------------------------------------
&namzgr_isf    !   isf cavity geometry definition                       (default: OFF)
!-----------------------------------------------------------------------
   rn_isfdep_min    = 10.         ! minimum isf draft tickness (if lower, isf draft set to this value)
   rn_glhw_min      = 1.e-3       ! minimum water column thickness to define the grounding line
   rn_isfhw_min     = 10          ! minimum water column thickness in the cavity once the grounding line defined.
   ln_isfchannel    = .false.     ! remove channel (based on 2d mask build from isfdraft-bathy)
   ln_isfconnect    = .false.     ! force connection under the ice shelf (based on 2d mask build from isfdraft-bathy)
      nn_kisfmax       = 999         ! limiter in level on the previous condition. (if change larger than this number, get back to value before we enforce the connection)
      rn_zisfmax       = 7000.       ! limiter in m     on the previous condition. (if change larger than this number, get back to value before we enforce the connection)
   ln_isfcheminey   = .false.     ! close cheminey
   ln_isfsubgl      = .false.     ! remove subglacial lake created by the remapping process
      rn_isfsubgllon   =    0.0      !  longitude of the seed to determine the open ocean
      rn_isfsubgllat   =    0.0      !  latitude  of the seed to determine the open ocean
/
  \end{forlines}
\end{listing}

  If the under ice shelf seas are opened (\np{ln_isfcav}{ln\_isfcav}), the depth of the ice shelf/ocean interface has to be included in 
  the \textit{isfdraft\_meter} file (Netcdf format). This file needs to include the \textit{isf\_draft} variable. 
  A positive value will mean ice shelf/ocean interface below the reference 0m ssh. 
  The exact shape of the ice shelf cavity (grounding line position and minimum thickness of the water column under an ice shelf, ...) can be specify in \nam{zgr_isf}{zgr\_isf} (\texttt{DOMAINcfg} only, \autoref{lst:namzgr_isf}). The details of each of these parameters is available in \autoref{subsec:ISF_dom} in the ISF section (\autoref{sec:LIO_isf}.
\section[Closed sea mask definition (\textit{domclo.F90})]{Closed sea mask definition (\protect\mdl{domclo})}
\begin{listing}
  \caption{\forcode{&namclo}}
  \label{lst:namdom_clo}
  \begin{forlines}
!-----------------------------------------------------------------------
&namclo ! (closed sea : need ln_domclo = .true. in namcfg)
!-----------------------------------------------------------------------
   rn_lon_opnsea = -2.0     ! longitude seed of open ocean
   rn_lat_opnsea = -2.0     ! latitude  seed of open ocean
   nn_closea = 8           ! number of closed seas ( = 0; only the open_sea mask will be computed)
   !                name   ! lon_src ! lat_src ! lon_trg ! lat_trg ! river mouth area   ! net evap/precip correction scheme ! radius tgt   ! id trg
   !                       ! (degree)! (degree)! (degree)! (degree)! local/coast/global ! (glo/rnf/emp)                     !     (m)      !
   ! North American lakes
   sn_lake(1) = 'superior' ,  -86.57 ,  47.30  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2
   sn_lake(2) = 'michigan' ,  -87.06 ,  42.74  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2
   sn_lake(3) = 'huron'    ,  -82.51 ,  44.74  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2
   sn_lake(4) = 'erie'     ,  -81.13 ,  42.25  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2
   sn_lake(5) = 'ontario'  ,  -77.72 ,  43.62  , -66.49  , 50.45   , 'local'            , 'rnf'                             ,   550000.0 , 2
   ! African Lake
   sn_lake(6) = 'victoria' ,   32.93 ,  -1.08  ,  30.44  , 31.37   , 'coast'            , 'emp'                             ,   100000.0 , 3
   ! Asian Lakes
   sn_lake(7) = 'caspian'  ,   50.0  ,  44.0   ,   0.0   ,  0.0    , 'global'           , 'glo'                             ,        0.0 , 1
   sn_lake(8) = 'aral'     ,   60.0  ,  45.0   ,   0.0   ,  0.0    , 'global'           , 'glo'                             ,        0.0 , 1
/
   \end{forlines}
\end{listing}

The options available to define the closed seas with \texttt{DOMAINcfg} are listed in \texttt{\&namclo}, while the control on how closed sea net fresh water input will be redistributed by NEMO is described in \autoref{subsec:DOM_closea} and \autoref{sec:MISC_closea}.
The individual definition of each closed sea is managed by \np{sn_lake}{sn\_lake}. In this fields the user needs to define:\\
   \begin{description}
   \item $\bullet$    The name of the closed sea (print output purposes).
   \item $\bullet$    The seed location to define the area of the closed sea (if seed on land because not present in this configuration, this closed sea will be ignored).
   \item $\bullet$    The seed location for the target area.
   \item $\bullet$    The type of target area ('local','coast' or 'global'). See point 6 for definition of these cases.
   \item $\bullet$    The type of redistribution scheme for the net fresh water flux over the closed sea (as a runoff in a target area, as emp in a target area, as emp globally). For the runoff case, if the net fwf is negative, it will be redistributed globally.
   \item $\bullet$    The radius of the target area (not used for the 'global' case). The target defined by a 'local' target area of a radius of 100km, for example, correspond sto all the wet points within this radius. The coastal case will return only the coastal point within the specified radius.
   \item $\bullet$    The target id. This id is used to group multiple lakes into the same river ouflow (Great Lakes for example).
   \end{description}

The closed sea module defines a number of masks in the \textit{domain\_cfg} output:
   \begin{description}
   \item[\textit{mask\_opensea}:] a mask of the main ocean without all the closed seas closed. This mask is defined by a flood filling algorithm with an initial seed (localisation defined by \np{rn_lon_opnsea}{rn\_lon\_opnsea} and \np{rn_lat_opnsea}{rn\_lat\_opnsea}).
   \item[\textit{mask\_csglo}, \textit{mask\_csrnf}, \textit{mask\_csemp}:] a mask of all the closed seas defined in the namelist by \np{sn_lake}{sn\_lake} for each redistribution scheme. The total number of defined closed seas has to be defined in \np{nn_closea}{nn\_closea}.
   \item[\textit{mask\_csgrpglo}, \textit{mask\_csgrprnf}, \textit{mask\_csgrpemp}:] a mask of all the closed seas and targets grouped by target id for each type of redistribution scheme.
   \item[\textit{mask\_csundef}:] a mask of all the closed sea not defined in \np{sn_lake}{sn\_lake}. This allows NEMO to mask them if needed or to inform the user of potential minor issues in its bathymetry.
   \end{description}

   
\subinc{\input{../../global/epilogue}}

\end{document}