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 11830 – NEMO

Changeset 11830


Ignore:
Timestamp:
2019-10-29T17:51:07+01:00 (4 years ago)
Author:
acc
Message:

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Merge changes to doc from trunk r10740 through r11740

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc
Files:
53 deleted
125 edited
65 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc

    • Property svn:ignore deleted
    • Property svn:mergeinfo deleted
    • Property svn:externals set to
      ^/utils/badges badges
      ^/utils/logos logos

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex

    • Property svn:ignore deleted

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/HTML_latex2html.sh

    r10146 r11830  
    11#!/bin/bash 
    22 
    3 ./inc/clean.sh 
    4 ./inc/build.sh 
     3#./inc/clean.sh 
     4#./inc/build.sh 
    55 
    6 cd tex_main 
    7 sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s/^/%/' \ 
    8    ../tex_sub/*.tex 
    9 sed -i    's#\\subfile{#\\include{#g' \ 
    10    NEMO_manual.tex 
    11 latex2html -local_icons -no_footnode -split 4 -link 2 -mkdir -dir ../html_LaTeX2HTML   \ 
    12             $*                                                                         \ 
    13    NEMO_manual 
    14 sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s/^%//' \ 
    15    ../tex_sub/*.tex 
    16 sed -i    's#\\include{#\\subfile{#g' \ 
    17    NEMO_manual.tex 
     6sed -i -e 's#utf8#latin1#'                                 \ 
     7       -e 's#\[outputdir=../build\]{minted}#\[\]{minted}#' \ 
     8       -e '/graphicspath/ s#{../#{../../#g'                \ 
     9          global/packages.tex 
     10 
     11cd ./NEMO/main 
     12sed -i -e 's#\\documentclass#%\\documentclass#' -e '/{document}/ s#^#%#'   ../subfiles/*.tex 
     13sed -i    's#\\subfile{#\\input{#'                                         chapters.tex appendices.tex 
     14 
     15#latex2html                                   -noimages -local_icons -no_footnode -split 4 -link 2 -dir ../html_LaTeX2HTML $* NEMO_manual 
     16latex2html -debug -noreuse -init_file ../../l2hconf.pm -local_icons                               -dir ../build/html               NEMO_manual 
     17 
     18sed -i -e 's#%\\documentclass#\\documentclass#' -e '/{document}/ s#^%##'   ../subfiles/*.tex 
     19sed -i    's#\\input{#\\subfile{#'                                         chapters.tex appendices.tex 
    1820cd - 
    1921 
     22sed -i -e 's#latin1#utf8#'                                 \ 
     23       -e 's#\[\]{minted}#\[outputdir=../build\]{minted}#' \ 
     24       -e '/graphicspath/ s#{../../#{../#g'                \ 
     25     global/packages.tex 
     26 
    2027exit 0 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO

    • Property svn:ignore deleted
    • Property svn:externals set to
      ^/utils/figures/NEMO figures

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/main

    • Property svn:ignore

      •  

        old new  
        1 *.aux 
        2 *.bbl 
        3 *.blg 
        4 *.dvi 
        5 *.fdb* 
        6 *.fls 
        7 *.idx 
        81*.ilg 
        92*.ind 
        10 *.log 
        11 *.maf 
        12 *.mtc* 
        13 *.out 
        14 *.pdf 
        15 *.toc 
        16 _minted-* 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/subfiles/chap_ASM.tex

    r10442 r11830  
    22 
    33\begin{document} 
    4 % ================================================================ 
    5 % Chapter Assimilation increments (ASM) 
    6 % ================================================================ 
     4 
    75\chapter{Apply Assimilation Increments (ASM)} 
    86\label{chap:ASM} 
    97 
    10 Authors: D. Lea,  M. Martin, K. Mogensen, A. Weaver, ...   % do we keep 
     8%    {\em 4.0} & {\em D. J. Lea} & {\em \NEMO\ 4.0 updates}  \\ 
     9%    {\em 3.4} & {\em D. J. Lea, M. Martin, K. Mogensen, A. Weaver} & {\em Initial version}  \\ 
    1110 
    12 \minitoc 
     11\thispagestyle{plain} 
    1312 
    14 \newpage 
     13\chaptertoc 
     14 
     15\paragraph{Changes record} ~\\ 
     16 
     17{\footnotesize 
     18  \begin{tabularx}{\textwidth}{l||X|X} 
     19    Release & Author(s) & Modifications \\ 
     20    \hline 
     21    {\em   4.0} & {\em ...} & {\em ...} \\ 
     22    {\em   3.6} & {\em ...} & {\em ...} \\ 
     23    {\em   3.4} & {\em ...} & {\em ...} \\ 
     24    {\em <=3.4} & {\em ...} & {\em ...} 
     25  \end{tabularx} 
     26} 
     27 
     28\clearpage 
    1529 
    1630The ASM code adds the functionality to apply increments to the model variables: temperature, salinity, 
    17 sea surface height, velocity and sea ice concentration.  
     31sea surface height, velocity and sea ice concentration. 
    1832These are read into the model from a NetCDF file which may be produced by separate data assimilation code. 
    1933The code can also output model background fields which are used as an input to data assimilation code. 
    20 This is all controlled by the namelist \textit{\ngn{nam\_asminc} }. 
     34This is all controlled by the namelist \nam{_asminc}{\_asminc}. 
    2135There is a brief description of all the namelist options provided. 
    2236To build the ASM code \key{asminc} must be set. 
    2337 
    24 %=============================================================== 
    25  
     38%% ================================================================================================= 
    2639\section{Direct initialization} 
    2740\label{sec:ASM_DI} 
    … …  
    2942Direct initialization (DI) refers to the instantaneous correction of the model background state using 
    3043the analysis increment. 
    31 DI is used when \np{ln\_asmdin} is set to true. 
     44DI is used when \np{ln_asmdin}{ln\_asmdin} is set to true. 
    3245 
     46%% ================================================================================================= 
    3347\section{Incremental analysis updates} 
    3448\label{sec:ASM_IAU} 
    … …  
    3751it may be preferable to introduce the increment gradually into the ocean model in order to 
    3852minimize spurious adjustment processes. 
    39 This technique is referred to as Incremental Analysis Updates (IAU) \citep{Bloom_al_MWR96}. 
     53This technique is referred to as Incremental Analysis Updates (IAU) \citep{bloom.takacs.ea_MWR96}. 
    4054IAU is a common technique used with 3D assimilation methods such as 3D-Var or OI. 
    41 IAU is used when \np{ln\_asmiau} is set to true. 
     55IAU is used when \np{ln_asmiau}{ln\_asmiau} is set to true. 
    4256 
    43 With IAU, the model state trajectory ${\bf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) 
     57With IAU, the model state trajectory ${\mathbf x}$ in the assimilation window ($t_{0} \leq t_{i} \leq t_{N}$) 
    4458is corrected by adding the analysis increments for temperature, salinity, horizontal velocity and SSH as 
    4559additional tendency terms to the prognostic equations: 
    4660\begin{align*} 
    47   % \label{eq:wa_traj_iau} 
    48   {\bf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\bf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\bf x}^{a} 
     61  % \label{eq:ASM_wa_traj_iau} 
     62  {\mathbf x}^{a}(t_{i}) = M(t_{i}, t_{0})[{\mathbf x}^{b}(t_{0})] \; + \; F_{i} \delta \tilde{\mathbf x}^{a} 
    4963\end{align*} 
    50 where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\bf x}^{a}$ defined such that 
     64where $F_{i}$ is a weighting function for applying the increments $\delta\tilde{\mathbf x}^{a}$ defined such that 
    5165$\sum_{i=1}^{N} F_{i}=1$. 
    52 ${\bf x}^b$ denotes the model initial state and ${\bf x}^a$ is the model state after the increments are applied. 
     66${\mathbf x}^b$ denotes the model initial state and ${\mathbf x}^a$ is the model state after the increments are applied. 
    5367To control the adjustment time of the model to the increment, 
    5468the increment can be applied over an arbitrary sub-window, $t_{m} \leq t_{i} \leq t_{n}$, 
    … …  
    5670Typically the increments are spread evenly over the full window. 
    5771In addition, two different weighting functions have been implemented. 
    58 The first function employs constant weights,  
     72The first function (namelist option \np{niaufn}{niaufn}=0) employs constant weights, 
    5973\begin{align} 
    60   \label{eq:F1_i} 
     74  \label{eq:ASM_F1_i} 
    6175  F^{(1)}_{i} 
    6276  =\left\{ 
    6377  \begin{array}{ll} 
    64     0     &    {\rm if} \; \; \; t_{i} < t_{m}                \\ 
    65     1/M &    {\rm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 
    66     0     &    {\rm if} \; \; \; t_{i} > t_{n} 
     78    0     &    {\mathrm if} \; \; \; t_{i} < t_{m}                \\ 
     79    1/M &    {\mathrm if} \; \; \; t_{m} < t_{i} \leq t_{n} \\ 
     80    0     &    {\mathrm if} \; \; \; t_{i} > t_{n} 
    6781  \end{array} 
    68             \right.  
     82            \right. 
    6983\end{align} 
    7084where $M = m-n$. 
    71 The second function employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, 
     85The second function (namelist option \np{niaufn}{niaufn}=1) employs peaked hat-like weights in order to give maximum weight in the centre of the sub-window, 
    7286with the weighting reduced linearly to a small value at the window end-points: 
    7387\begin{align} 
    74   \label{eq:F2_i} 
     88  \label{eq:ASM_F2_i} 
    7589  F^{(2)}_{i} 
    7690  =\left\{ 
    7791  \begin{array}{ll} 
    78     0                           &    {\rm if} \; \; \; t_{i}       <     t_{m}                        \\ 
    79     \alpha \, i               &    {\rm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ 
    80     \alpha \, (M - i +1) &    {\rm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 
    81     0                            &   {\rm if} \; \; \; t_{i}        >    t_{n} 
     92    0                           &    {\mathrm if} \; \; \; t_{i}       <     t_{m}                        \\ 
     93    \alpha \, i               &    {\mathrm if} \; \; \; t_{m}    \leq t_{i}    \leq   t_{M/2}   \\ 
     94    \alpha \, (M - i +1) &    {\mathrm if} \; \; \; t_{M/2}  <    t_{i}    \leq   t_{n}       \\ 
     95    0                            &   {\mathrm if} \; \; \; t_{i}        >    t_{n} 
    8296  \end{array} 
    8397                                   \right. 
    8498\end{align} 
    85 where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even.  
    86 The weights described by \autoref{eq:F2_i} provide a smoother transition of the analysis trajectory from 
    87 one assimilation cycle to the next than that described by \autoref{eq:F1_i}. 
     99where $\alpha^{-1} = \sum_{i=1}^{M/2} 2i$ and $M$ is assumed to be even. 
     100The weights described by \autoref{eq:ASM_F2_i} provide a smoother transition of the analysis trajectory from 
     101one assimilation cycle to the next than that described by \autoref{eq:ASM_F1_i}. 
    88102 
    89 %========================================================================== 
    90 % Divergence damping description %%% 
     103%% ================================================================================================= 
    91104\section{Divergence damping initialisation} 
    92105\label{sec:ASM_div_dmp} 
    93106 
    94 The velocity increments may be initialized by the iterative application of a divergence damping operator. 
    95 In iteration step $n$ new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 
     107It is quite challenging for data assimilation systems to provide non-divergent velocity increments. 
     108Applying divergent velocity increments will likely cause spurious vertical velocities in the model. This section describes a method to take velocity increments provided to \NEMO\ ($u^0_I$ and $v^0_I$) and adjust them by the iterative application of a divergence damping operator. The method is also described in \citet{dobricic.pinardi.ea_OS07}. 
     109 
     110In iteration step $n$ (starting at $n=1$) new estimates of velocity increments $u^{n}_I$ and $v^{n}_I$ are updated by: 
     111 
    96112\begin{equation} 
    97   \label{eq:asm_dmp} 
     113  \label{eq:ASM_dmp} 
    98114  \left\{ 
    99115    \begin{aligned} 
    … …  
    105121  \right., 
    106122\end{equation} 
    107 where 
     123 
     124where the divergence is defined as 
     125 
    108126\[ 
    109   % \label{eq:asm_div} 
     127  % \label{eq:ASM_div} 
    110128  \chi^{n-1}_I = \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    111129  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u^{n-1}_I} \right] 
    112130      +\delta_j \left[ {e_{1v}\,e_{3v}\,v^{n-1}_I} \right]} \right). 
    113131\] 
    114 By the application of \autoref{eq:asm_dmp} and \autoref{eq:asm_dmp} the divergence is filtered in each iteration, 
     132 
     133By the application of \autoref{eq:ASM_dmp} the divergence is filtered in each iteration, 
    115134and the vorticity is left unchanged. 
    116135In the presence of coastal boundaries with zero velocity increments perpendicular to the coast 
    … …  
    118137This type of the initialisation reduces the vertical velocity magnitude and 
    119138alleviates the problem of the excessive unphysical vertical mixing in the first steps of the model integration 
    120 \citep{Talagrand_JAS72, Dobricic_al_OS07}. 
     139\citep{talagrand_JAS72, dobricic.pinardi.ea_OS07}. 
    121140Diffusion coefficients are defined as $A_D = \alpha e_{1t} e_{2t}$, where $\alpha = 0.2$. 
    122 The divergence damping is activated by assigning to \np{nn\_divdmp} in the \textit{nam\_asminc} namelist 
     141The divergence damping is activated by assigning to \np{nn_divdmp}{nn\_divdmp} in the \nam{_asminc}{\_asminc} namelist 
    123142a value greater than zero. 
    124 By choosing this value to be of the order of 100 the increments in 
    125 the vertical velocity will be significantly reduced. 
     143This specifies the number of iterations of the divergence damping. Setting a value of the order of 100 will result in a significant reduction in the vertical velocity induced by the increments. 
    126144 
    127  
    128 %========================================================================== 
    129  
     145%% ================================================================================================= 
    130146\section{Implementation details} 
    131147\label{sec:ASM_details} 
    132148 
    133 Here we show an example \ngn{namasm} namelist and the header of an example assimilation increments file on 
     149Here we show an example \nam{_asminc}{\_asminc} namelist and the header of an example assimilation increments file on 
    134150the ORCA2 grid. 
    135151 
    136 %------------------------------------------namasm----------------------------------------------------- 
    137 % 
    138 \nlst{nam_asminc} 
    139 %------------------------------------------------------------------------------------------------------------- 
     152\begin{listing} 
     153  \nlst{nam_asminc} 
     154  \caption{\forcode{&nam_asminc}} 
     155  \label{lst:nam_asminc} 
     156\end{listing} 
    140157 
    141158The header of an assimilation increments file produced using the NetCDF tool 
    … …  
    177194\end{clines} 
    178195 
    179 \biblio 
    180  
    181 \pindex 
     196\subinc{\input{../../global/epilogue}} 
    182197 
    183198\end{document} 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/subfiles/chap_DIA.tex

    r10509 r11830  
    22 
    33\begin{document} 
    4 % ================================================================ 
    5 % Chapter I/O & Diagnostics 
    6 % ================================================================ 
     4 
    75\chapter{Output and Diagnostics (IOM, DIA, TRD, FLO)} 
    86\label{chap:DIA} 
    97 
    10 \minitoc 
    11  
    12 \newpage 
    13  
    14 % ================================================================ 
    15 %       Old Model Output  
    16 % ================================================================ 
    17 \section{Old model output (default)} 
     8%    {\em 4.0} & {\em Mirek Andrejczuk, Massimiliano Drudi} & {\em }  \\ 
     9%    {\em }      & {\em Dorotea Iovino, Nicolas Martin} & {\em }  \\ 
     10%    {\em 3.6} & {\em Gurvan Madec, Sebastien Masson } & {\em }  \\ 
     11%    {\em 3.4} & {\em Gurvan Madec, Rachid Benshila, Andrew Coward } & {\em }  \\ 
     12%    {\em }      & {\em Christian Ethe, Sebastien Masson } & {\em }  \\ 
     13 
     14\thispagestyle{plain} 
     15 
     16\chaptertoc 
     17 
     18\paragraph{Changes record} ~\\ 
     19 
     20{\footnotesize 
     21  \begin{tabularx}{\textwidth}{l||X|X} 
     22    Release & Author(s) & Modifications \\ 
     23    \hline 
     24    {\em   4.0} & {\em ...} & {\em ...} \\ 
     25    {\em   3.6} & {\em ...} & {\em ...} \\ 
     26    {\em   3.4} & {\em ...} & {\em ...} \\ 
     27    {\em <=3.4} & {\em ...} & {\em ...} 
     28  \end{tabularx} 
     29} 
     30 
     31\clearpage 
     32 
     33%% ================================================================================================= 
     34\section{Model output} 
    1835\label{sec:DIA_io_old} 
    1936 
    … …  
    2542the same run performed in one step. 
    2643It should be noted that this requires that the restart file contains two consecutive time steps for 
    27 all the prognostic variables, and that it is saved in the same binary format as the one used by the computer that 
    28 is to read it (in particular, 32 bits binary IEEE format must not be used for this file). 
     44all the prognostic variables. 
    2945 
    3046The output listing and file(s) are predefined but should be checked and eventually adapted to the user's needs. 
    31 The output listing is stored in the $ocean.output$ file. 
    32 The information is printed from within the code on the logical unit $numout$. 
     47The output listing is stored in the \textit{ocean.output} file. 
     48The information is printed from within the code on the logical unit \texttt{numout}. 
    3349To locate these prints, use the UNIX command "\textit{grep -i numout}" in the source code directory. 
    3450 
    … …  
    3955A complete description of the use of this I/O server is presented in the next section. 
    4056 
    41 By default, \key{iomput} is not defined, 
    42 NEMO produces NetCDF with the old IOIPSL library which has been kept for compatibility and its easy installation. 
    43 However, the IOIPSL library is quite inefficient on parallel machines and, since version 3.2, 
    44 many diagnostic options have been added presuming the use of \key{iomput}. 
    45 The usefulness of the default IOIPSL-based option is expected to reduce with each new release. 
    46 If \key{iomput} is not defined, output files and content are defined in the \mdl{diawri} module and 
    47 contain mean (or instantaneous if \key{diainstant} is defined) values over a regular period of 
    48 nn\_write time-steps (namelist parameter).  
    49  
    50 %\gmcomment{                    % start of gmcomment 
    51  
    52 % ================================================================ 
    53 % Diagnostics 
    54 % ================================================================ 
     57%\cmtgm{                    % start of gmcomment 
     58 
     59%% ================================================================================================= 
    5560\section{Standard model output (IOM)} 
    5661\label{sec:DIA_iom} 
    5762 
    58 Since version 3.2, iomput is the NEMO output interface of choice. 
     63Since version 3.2, iomput is the \NEMO\ output interface of choice. 
    5964It has been designed to be simple to use, flexible and efficient. 
    60 The two main purposes of iomput are:  
     65The two main purposes of iomput are: 
    6166 
    6267\begin{enumerate} 
    63 \item 
    64   The complete and flexible control of the output files through external XML files adapted by 
     68\item The complete and flexible control of the output files through external XML files adapted by 
    6569  the user from standard templates. 
    66 \item 
    67   To achieve high performance and scalable output through the optional distribution of 
     70\item To achieve high performance and scalable output through the optional distribution of 
    6871  all diagnostic output related tasks to dedicated processes. 
    6972\end{enumerate} 
    7073 
    71 The first functionality allows the user to specify, without code changes or recompilation,  
     74The first functionality allows the user to specify, without code changes or recompilation, 
    7275aspects of the diagnostic output stream, such as: 
    7376 
    7477\begin{itemize} 
    75 \item 
    76   The choice of output frequencies that can be different for each file (including real months and years). 
    77 \item 
    78   The choice of file contents; includes complete flexibility over which data are written in which files 
     78\item The choice of output frequencies that can be different for each file (including real months and years). 
     79\item The choice of file contents; includes complete flexibility over which data are written in which files 
    7980  (the same data can be written in different files). 
    80 \item 
    81   The possibility to split output files at a chosen frequency. 
    82 \item 
    83   The possibility to extract a vertical or an horizontal subdomain. 
    84 \item 
    85   The choice of the temporal operation to perform, \eg: average, accumulate, instantaneous, min, max and once. 
    86 \item 
    87   Control over metadata via a large XML "database" of possible output fields. 
     81\item The possibility to split output files at a chosen frequency. 
     82\item The possibility to extract a vertical or an horizontal subdomain. 
     83\item The choice of the temporal operation to perform, \eg: average, accumulate, instantaneous, min, max and once. 
     84\item Control over metadata via a large XML "database" of possible output fields. 
    8885\end{itemize} 
    8986 
    … …  
    9188in a very easy way. 
    9289All details of iomput functionalities are listed in the following subsections. 
    93 Examples of the XML files that control the outputs can be found in: \path{NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef.xml}, 
    94 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml}. \\ 
     90Examples of the XML files that control the outputs can be found in: 
     91\path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef.xml}, 
     92\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     93\path{cfgs/SHARED/field_def_nemo-pisces.xml}, 
     94\path{cfgs/SHARED/field_def_nemo-ice.xml} and \path{cfgs/SHARED/domain_def_nemo.xml}. \\ 
    9595 
    9696The second functionality targets output performance when running in parallel (\key{mpp\_mpi}). 
    97 Iomput provides the possibility to specify N dedicated I/O processes (in addition to the NEMO processes) 
     97Iomput provides the possibility to specify N dedicated I/O processes (in addition to the \NEMO\ processes) 
    9898to collect and write the outputs. 
    9999With an appropriate choice of N by the user, the bottleneck associated with the writing of 
    100100the output files can be greatly reduced. 
    101101 
    102 In version 3.6, the iom\_put interface depends on 
    103 an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-1.0}{XIOS-1.0}  
    104 (use of revision 618 or higher is required). 
     102In version 3.6, the \rou{iom\_put} interface depends on 
     103an external code called \href{https://forge.ipsl.jussieu.fr/ioserver/browser/XIOS/branchs/xios-2.5}{XIOS-2.5} 
     104%(use of revision 618 or higher is required). 
    105105This new IO server can take advantage of the parallel I/O functionality of NetCDF4 to 
    106106create a single output file and therefore to bypass the rebuilding phase. 
    107107Note that writing in parallel into the same NetCDF files requires that your NetCDF4 library is linked to 
    108 an HDF5 library that has been correctly compiled (\ie with the configure option $--$enable-parallel). 
     108an HDF5 library that has been correctly compiled (\ie\ with the configure option $--$enable-parallel). 
    109109Note that the files created by iomput through XIOS are incompatible with NetCDF3. 
    110110All post-processsing and visualization tools must therefore be compatible with NetCDF4 and not only NetCDF3. 
    111111 
    112112Even if not using the parallel I/O functionality of NetCDF4, using N dedicated I/O servers, 
    113 where N is typically much less than the number of NEMO processors, will reduce the number of output files created. 
    114 This can greatly reduce the post-processing burden usually associated with using large numbers of NEMO processors. 
     113where N is typically much less than the number of \NEMO\ processors, will reduce the number of output files created. 
     114This can greatly reduce the post-processing burden usually associated with using large numbers of \NEMO\ processors. 
    115115Note that for smaller configurations, the rebuilding phase can be avoided, 
    116116even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server. 
    117117 
     118%% ================================================================================================= 
    118119\subsection{XIOS: Reading and writing restart file} 
    119120 
    120 XIOS may be used to read single file restart produced by NEMO. Currently only the variables written to  
    121 file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np{ln\_xios\_read}\forcode{ = .true. } 
    122 in \textit{namelist\_cfg}. This setting will be ignored when multiple restart files are present, and default NEMO  
    123 functionality will be used for reading. There is no need to change iodef.xml file to use XIOS to read  
    124 restart, all definitions are done within the NEMO code. For high resolution configuration, however,  
     121XIOS may be used to read single file restart produced by \NEMO. Currently only the variables written to 
     122file \forcode{numror} can be handled by XIOS. To activate restart reading using XIOS, set \np[=.true. ]{ln_xios_read}{ln\_xios\_read} 
     123in \textit{namelist\_cfg}. This setting will be ignored when multiple restart files are present, and default \NEMO 
     124functionality will be used for reading. There is no need to change iodef.xml file to use XIOS to read 
     125restart, all definitions are done within the \NEMO\ code. For high resolution configuration, however, 
    125126there may be a need to add the following line in iodef.xml (xios context): 
    126127 
    … …  
    129130\end{xmllines} 
    130131 
    131 This variable sets timeout for reading.  
    132  
    133 If XIOS is to be used to read restart from file generated with an earlier NEMO version (3.6 for instance), 
     132This variable sets timeout for reading. 
     133 
     134If XIOS is to be used to read restart from file generated with an earlier \NEMO\ version (3.6 for instance), 
    134135dimension \forcode{z} defined in restart file must be renamed to \forcode{nav_lev}.\\ 
    135136 
    136 XIOS can also be used to write NEMO restart. A namelist parameter \np{nn\_wxios} is used to determine the  
    137 type of restart NEMO will write. If it is set to 0, default NEMO functionality will be used - each  
    138 processor writes its own restart file; if it is set to 1 XIOS will write restart into a single file;  
    139 for \np{nn\_wxios = 2} the restart will be written by XIOS into multiple files, one for each XIOS server.  
    140 Note, however, that \textbf{NEMO will not read restart generated by XIOS when \np{nn\_wxios = 2}}. The restart will  
    141 have to be rebuild before continuing the run. This option aims to reduce number of restart files generated by NEMO only,  
    142 and may be useful when there is a need to change number of processors used to run simulation.  
     137XIOS can also be used to write \NEMO\ restart. A namelist parameter \np{nn_wxios}{nn\_wxios} is used to determine the 
     138type of restart \NEMO\ will write. If it is set to 0, default \NEMO\ functionality will be used - each 
     139processor writes its own restart file; if it is set to 1 XIOS will write restart into a single file; 
     140for \np[=2]{nn_wxios}{nn\_wxios} the restart will be written by XIOS into multiple files, one for each XIOS server. 
     141Note, however, that \textbf{\NEMO\ will not read restart generated by XIOS when \np[=2]{nn_wxios}{nn\_wxios}}. The restart will 
     142have to be rebuild before continuing the run. This option aims to reduce number of restart files generated by \NEMO\ only, 
     143and may be useful when there is a need to change number of processors used to run simulation. 
    143144 
    144145If an additional variable must be written to a restart file, the following steps are needed: 
    145 \begin{description} 
    146    \item[step 1:] add variable name to a list of restart variables (in subroutine \rou{iom\_set\_rst\_vars,} \mdl{iom}) and  
    147 define correct grid for the variable (\forcode{grid_N_3D} - 3D variable, \forcode{grid_N} - 2D variable, \forcode{grid_vector} -  
     146\begin{enumerate} 
     147\item Add variable name to a list of restart variables (in subroutine \rou{iom\_set\_rst\_vars,} \mdl{iom}) and 
     148define correct grid for the variable (\forcode{grid_N_3D} - 3D variable, \forcode{grid_N} - 2D variable, \forcode{grid_vector} - 
    1481491D variable, \forcode{grid_scalar} - scalar), 
    149    \item[step 2:] add variable to the list of fields written by restart.  This can be done either in subroutine  
    150 \rou{iom\_set\_rstw\_core} (\mdl{iom}) or by calling  \rou{iom\_set\_rstw\_active} (\mdl{iom}) with the name of a variable  
    151 as an argument. This convention follows approach for writing restart using iom, where variables are  
     150\item Add variable to the list of fields written by restart.  This can be done either in subroutine 
     151\rou{iom\_set\_rstw\_core} (\mdl{iom}) or by calling  \rou{iom\_set\_rstw\_active} (\mdl{iom}) with the name of a variable 
     152as an argument. This convention follows approach for writing restart using iom, where variables are 
    152153written either by \rou{rst\_write} or by calling \rou{iom\_rstput} from individual routines. 
    153 \end{description} 
    154  
     154\end{enumerate} 
    155155 
    156156An older versions of XIOS do not support reading functionality. It's recommended to use at least XIOS2@1451. 
    157157 
    158  
     158%% ================================================================================================= 
    159159\subsection{XIOS: XML Inputs-Outputs Server} 
    160160 
     161%% ================================================================================================= 
    161162\subsubsection{Attached or detached mode?} 
    162163 
    … …  
    168169\xmlline|<variable id="using_server" type="bool"></variable>| 
    169170 
    170 The {\tt using\_server} setting determines whether or not the server will be used in \textit{attached mode} 
    171 (as a library) [{\tt> false <}] or in \textit{detached mode} 
    172 (as an external executable on N additional, dedicated cpus) [{\tt > true <}]. 
    173 The \textit{attached mode} is simpler to use but much less efficient for massively parallel applications. 
     171The \texttt{using\_server} setting determines whether or not the server will be used in 
     172\textit{attached mode} 
     173(as a library) [\texttt{> false <}] or in \textit{detached mode} 
     174(as an external executable on N additional, dedicated cpus) [\texttt{ > true <}]. 
     175The \textit{attached mode} is simpler to use but much less efficient for 
     176massively parallel applications. 
    174177The type of each file can be either ''multiple\_file'' or ''one\_file''. 
    175178 
    176179In \textit{attached mode} and if the type of file is ''multiple\_file'', 
    177 then each NEMO process will also act as an IO server and produce its own set of output files. 
     180then each \NEMO\ process will also act as an IO server and produce its own set of output files. 
    178181Superficially, this emulates the standard behaviour in previous versions. 
    179182However, the subdomain written out by each process does not correspond to 
    … …  
    187190write to its own set of output files. 
    188191If the ''one\_file'' option is chosen then all XIOS processes will collect their longitudinal strips and 
    189 write (in parallel) to a single output file.  
     192write (in parallel) to a single output file. 
    190193Note running in detached mode requires launching a Multiple Process Multiple Data (MPMD) parallel job. 
    191194The following subsection provides a typical example but the syntax will vary in different MPP environments. 
    192195 
     196%% ================================================================================================= 
    193197\subsubsection{Number of cpu used by XIOS in detached mode} 
    194198 
    195199The number of cores used by the XIOS is specified when launching the model. 
    196 The number of cores dedicated to XIOS should be from \texttildelow1/10 to \texttildelow1/50 of the number of  
    197 cores dedicated to NEMO. 
     200The number of cores dedicated to XIOS should be from \texttildelow1/10 to \texttildelow1/50 of the number of 
     201cores dedicated to \NEMO. 
    198202Some manufacturers suggest using O($\sqrt{N}$) dedicated IO processors for N processors but 
    199 this is a general recommendation and not specific to NEMO. 
     203this is a general recommendation and not specific to \NEMO. 
    200204It is difficult to provide precise recommendations because the optimal choice will depend on 
    201 the particular hardware properties of the target system  
     205the particular hardware properties of the target system 
    202206(parallel filesystem performance, available memory, memory bandwidth etc.) 
    203207and the volume and frequency of data to be created. 
    … …  
    205209\cmd|mpirun -np 62 ./nemo.exe : -np 2 ./xios_server.exe| 
    206210 
     211%% ================================================================================================= 
    207212\subsubsection{Control of XIOS: the context in iodef.xml} 
    208213 
    209 As well as the {\tt using\_server} flag, other controls on the use of XIOS are set in the XIOS context in iodef.xml. 
     214As well as the \texttt{using\_server} flag, other controls on the use of XIOS are set in 
     215the XIOS context in \textit{iodef.xml}. 
    210216See the XML basics section below for more details on XML syntax and rules. 
    211217 
    212218\begin{table} 
    213   \scriptsize 
    214219  \begin{tabularx}{\textwidth}{|lXl|} 
    215220    \hline 
    … …  
    220225    \hline 
    221226    buffer\_size                                                            & 
    222     buffer size used by XIOS to send data from NEMO to XIOS. 
     227    buffer size used by XIOS to send data from \NEMO\ to XIOS. 
    223228    Larger is more efficient. 
    224229    Note that needed/used buffer sizes are summarized at the end of the job & 
    … …  
    226231    \hline 
    227232    buffer\_server\_factor\_size                                            & 
    228     ratio between NEMO and XIOS buffer size. 
     233    ratio between \NEMO\ and XIOS buffer size. 
    229234    Should be 2.                                                            & 
    230235    2        \\ 
    … …  
    243248    \hline 
    244249    oasis\_codes\_id                                                        & 
    245     when using oasis, define the identifier of NEMO in the namcouple. 
     250    when using oasis, define the identifier of \NEMO\ in the namcouple. 
    246251    Note that the identifier of XIOS is xios.x                              & 
    247252    oceanx   \\ 
    … …  
    250255\end{table} 
    251256 
     257%% ================================================================================================= 
    252258\subsection{Practical issues} 
    253259 
     260%% ================================================================================================= 
    254261\subsubsection{Installation} 
    255262 
    256 As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with NEMO. 
     263As mentioned, XIOS is supported separately and must be downloaded and compiled before it can be used with \NEMO. 
    257264See the installation guide on the \href{http://forge.ipsl.jussieu.fr/ioserver/wiki}{XIOS} wiki for help and guidance. 
    258 NEMO will need to link to the compiled XIOS library. 
    259 The \href{https://forge.ipsl.jussieu.fr/nemo/wiki/Users/ModelInterfacing/InputsOutputs#Inputs-OutputsusingXIOS} 
    260 {XIOS with NEMO} guide provides an example illustration of how this can be achieved. 
    261  
     265\NEMO\ will need to link to the compiled XIOS library. 
     266The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/install.html#extract-and-install-xios} 
     267{Extract and install XIOS} guide provides an example illustration of how this can be achieved. 
     268 
     269%% ================================================================================================= 
    262270\subsubsection{Add your own outputs} 
    263271 
    … …  
    268276 
    269277\begin{enumerate} 
    270 \item[1.] 
    271   in NEMO code, add a \forcode{CALL iom\_put( 'identifier', array )} where you want to output a 2D or 3D array. 
    272 \item[2.] 
    273   If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in 
     278\item in \NEMO\ code, add a \forcode{CALL iom_put( 'identifier', array )} where you want to output a 2D or 3D array. 
     279\item If necessary, add \forcode{USE iom ! I/O manager library} to the list of used modules in 
    274280  the upper part of your module. 
    275 \item[3.] 
    276   in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code 
     281\item in the field\_def.xml file, add the definition of your variable using the same identifier you used in the f90 code 
    277282  (see subsequent sections for a details of the XML syntax and rules). 
    278283  For example: 
    279  
    280284\begin{xmllines} 
    281285<field_definition> 
    282286   <field_group id="grid_T" grid_ref="grid_T_3D"> <!-- T grid --> 
    283287   ... 
    284       <field id="identifier" long_name="blabla" ... />    
     288      <field id="identifier" long_name="blabla" ... /> 
    285289      ... 
    286 </field_definition>  
    287 \end{xmllines} 
    288  
     290</field_definition> 
     291\end{xmllines} 
    289292Note your definition must be added to the field\_group whose reference grid is consistent with the size of 
    290293the array passed to iomput. 
    … …  
    293296(iom\_set\_domain\_attr and iom\_set\_axis\_attr in \mdl{iom}) or defined in the domain\_def.xml file. 
    294297\eg: 
    295  
    296298\begin{xmllines} 
    297299<grid id="grid_T_3D" domain_ref="grid_T" axis_ref="deptht"/> 
    298300\end{xmllines} 
    299  
    300 Note, if your array is computed within the surface module each \np{nn\_fsbc} time\_step, 
     301Note, if your array is computed within the surface module each \np{nn_fsbc}{nn\_fsbc} time\_step, 
    301302add the field definition within the field\_group defined with the id "SBC": 
    302303\xmlcode{<field_group id="SBC" ...>} which has been defined with the correct frequency of operations 
    303304(iom\_set\_field\_attr in \mdl{iom}) 
    304 \item[4.] 
    305   add your field in one of the output files defined in iodef.xml 
     305\item add your field in one of the output files defined in iodef.xml 
    306306  (again see subsequent sections for syntax and rules) 
    307  
    308 \begin{xmllines} 
    309 <file id="file1" .../>    
     307\begin{xmllines} 
     308<file id="file1" .../> 
    310309... 
    311    <field field_ref="identifier" />    
     310   <field field_ref="identifier" /> 
    312311   ... 
    313 </file>    
    314 \end{xmllines} 
    315  
     312</file> 
     313\end{xmllines} 
    316314\end{enumerate} 
    317315 
     316%% ================================================================================================= 
    318317\subsection{XML fundamentals} 
    319318 
     319%% ================================================================================================= 
    320320\subsubsection{ XML basic rules} 
    321321 
    … …  
    327327See \href{http://www.xmlnews.org/docs/xml-basics.html}{here} for more details. 
    328328 
    329 \subsubsection{Structure of the XML file used in NEMO} 
     329%% ================================================================================================= 
     330\subsubsection{Structure of the XML file used in \NEMO} 
    330331 
    331332The XML file used in XIOS is structured by 7 families of tags: 
    … …  
    334335 
    335336\begin{table} 
    336   \scriptsize 
    337337  \begin{tabular*}{\textwidth}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} 
    338338    \hline 
    … …  
    366366 
    367367\begin{table} 
    368   \scriptsize 
    369368  \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} 
    370369    \hline 
    … …  
    376375                                                                                                     \xmlcode{<context id="xios" ... >}   \\ 
    377376    \hline 
    378     context nemo    &   context containing IO information for NEMO (mother grid when using AGRIF)  &  
     377    context nemo    &   context containing IO information for \NEMO\ (mother grid when using AGRIF)  & 
    379378                                                                                                     \xmlcode{<context id="nemo" ... >}   \\ 
    380379    \hline 
    381     context 1\_nemo &   context containing IO information for NEMO child grid 1 (when using AGRIF) &  
     380    context 1\_nemo &   context containing IO information for \NEMO\ child grid 1 (when using AGRIF) & 
    382381                                                                                                     \xmlcode{<context id="1_nemo" ... >} \\ 
    383382    \hline 
    384     context n\_nemo &   context containing IO information for NEMO child grid n (when using AGRIF) &  
     383    context n\_nemo &   context containing IO information for \NEMO\ child grid n (when using AGRIF) & 
    385384                                                                                                     \xmlcode{<context id="n_nemo" ... >} \\ 
    386385    \hline 
    … …  
    391390 
    392391\begin{table} 
    393   \scriptsize 
    394392  \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} 
    395393    \hline 
    … …  
    407405\end{table} 
    408406 
    409 \noindent Each context tag related to NEMO (mother or child grids) is divided into 5 parts  
     407\noindent Each context tag related to \NEMO\ (mother or child grids) is divided into 5 parts 
    410408(that can be defined in any order): 
    411409 
    412410\begin{table} 
    413   \scriptsize 
    414411  \begin{tabular}{|p{0.15\textwidth}p{0.4\textwidth}p{0.35\textwidth}|} 
    415412    \hline 
    … …  
    436433\end{table} 
    437434 
     435%% ================================================================================================= 
    438436\subsubsection{Nesting XML files} 
    439437 
    … …  
    441439The inclusion of XML files into the main XML file can be done through the attribute src: 
    442440\xmlline|<context src="./nemo_def.xml" />| 
    443   
    444 \noindent In NEMO, by default, the field and domain definition is done in 2 separate files: 
    445 \path{NEMOGCM/CONFIG/SHARED/field_def.xml} and \path{NEMOGCM/CONFIG/SHARED/domain_def.xml} that 
     441 
     442\noindent In \NEMO, by default, the field definition is done in 3 separate files ( 
     443\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     444\path{cfgs/SHARED/field_def_nemo-pisces.xml} and 
     445\path{cfgs/SHARED/field_def_nemo-ice.xml} ) and the  domain definition is done in another file ( \path{cfgs/SHARED/domain_def_nemo.xml} ) 
     446that 
    446447are included in the main iodef.xml file through the following commands: 
    447448\begin{xmllines} 
    448 <field_definition src="./field_def.xml" /> 
    449 <domain_definition src="./domain_def.xml" /> 
    450 \end{xmllines} 
    451  
     449<context id="nemo" src="./context_nemo.xml"/> 
     450\end{xmllines} 
     451 
     452%% ================================================================================================= 
    452453\subsubsection{Use of inheritance} 
    453454 
    … …  
    462463\begin{xmllines} 
    463464<field_definition operation="average" > 
    464    <field id="sst"                    />   <!-- averaged      sst -->  
    465    <field id="sss" operation="instant"/>   <!-- instantaneous sss -->  
    466 </field_definition>  
     465   <field id="sst"                    />   <!-- averaged      sst --> 
     466   <field id="sss" operation="instant"/>   <!-- instantaneous sss --> 
     467</field_definition> 
    467468\end{xmllines} 
    468469 
    … …  
    481482</field_definition> 
    482483<file_definition> 
    483    <file id="myfile" output_freq="1d" />    
     484   <file id="myfile" output_freq="1d" /> 
    484485      <field field_ref="sst"                            />  <!-- default def --> 
    485486      <field field_ref="sss" long_name="my description" />  <!-- overwrite   --> 
    486487   </file> 
    487 </file_definition>  
     488</file_definition> 
    488489\end{xmllines} 
    489490 
    490491Inherit (and overwrite, if needed) the attributes of a tag you are refering to. 
    491492 
     493%% ================================================================================================= 
    492494\subsubsection{Use of groups} 
    493495 
    … …  
    508510 
    509511Secondly, the group can be used to replace a list of elements. 
    510 Several examples of groups of fields are proposed at the end of the file \path{CONFIG/SHARED/field_def.xml}. 
     512Several examples of groups of fields are proposed at the end of the XML field files ( 
     513\path{cfgs/SHARED/field_def_nemo-oce.xml}, 
     514\path{cfgs/SHARED/field_def_nemo-pisces.xml} and 
     515\path{cfgs/SHARED/field_def_nemo-ice.xml} ) . 
    511516For example, a short list of the usual variables related to the U grid: 
    512517 
    … …  
    514519<field_group id="groupU" > 
    515520   <field field_ref="uoce"  /> 
    516    <field field_ref="suoce" /> 
     521   <field field_ref="ssu" /> 
    517522   <field field_ref="utau"  /> 
    518523</field_group> 
    … …  
    525530   <field_group group_ref="groupU" /> 
    526531   <field field_ref="uocetr_eff"   />  <!-- add another field --> 
    527 </file>    
    528 \end{xmllines} 
    529  
     532</file> 
     533\end{xmllines} 
     534 
     535%% ================================================================================================= 
    530536\subsection{Detailed functionalities} 
    531537 
    … …  
    533539the new functionalities offered by the XML interface of XIOS. 
    534540 
     541%% ================================================================================================= 
    535542\subsubsection{Define horizontal subdomains} 
    536543 
    537544Horizontal subdomains are defined through the attributs zoom\_ibegin, zoom\_jbegin, zoom\_ni, zoom\_nj of 
    538545the tag family domain. 
    539 It must therefore be done in the domain part of the XML file.  
    540 For example, in \path{CONFIG/SHARED/domain_def.xml}, we provide the following example of a definition of  
     546It must therefore be done in the domain part of the XML file. 
     547For example, in \path{cfgs/SHARED/domain_def.xml}, we provide the following example of a definition of 
    541548a 5 by 5 box with the bottom left corner at point (10,10). 
    542549 
    543550\begin{xmllines} 
    544 <domain_group id="grid_T"> 
    545    <domain id="myzoom" zoom_ibegin="10" zoom_jbegin="10" zoom_ni="5" zoom_nj="5" /> 
     551<domain id="myzoomT" domain_ref="grid_T"> 
     552   <zoom_domain ibegin="10" jbegin="10" ni="5" nj="5" /> 
    546553\end{xmllines} 
    547554 
    … …  
    551558\begin{xmllines} 
    552559<file id="myfile_vzoom" output_freq="1d" > 
    553    <field field_ref="toce" domain_ref="myzoom"/> 
     560   <field field_ref="toce" domain_ref="myzoomT"/> 
    554561</file> 
    555562\end{xmllines} 
    556563 
    557 Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain.  
     564Moorings are seen as an extrem case corresponding to a 1 by 1 subdomain. 
    558565The Equatorial section, the TAO, RAMA and PIRATA moorings are already registered in the code and 
    559566can therefore be outputted without taking care of their (i,j) position in the grid. 
    … …  
    568575\end{xmllines} 
    569576 
    570 Note that if the domain decomposition used in XIOS cuts the subdomain in several parts and if  
    571 you use the ''multiple\_file'' type for your output files,  
    572 you will endup with several files you will need to rebuild using unprovided tools (like ncpdq and ncrcat,  
     577Note that if the domain decomposition used in XIOS cuts the subdomain in several parts and if 
     578you use the ''multiple\_file'' type for your output files, 
     579you will endup with several files you will need to rebuild using unprovided tools (like ncpdq and ncrcat, 
    573580\href{http://nco.sourceforge.net/nco.html#Concatenation}{see nco manual}). 
    574581We are therefore advising to use the ''one\_file'' type in this case. 
    575582 
     583%% ================================================================================================= 
    576584\subsubsection{Define vertical zooms} 
    577585 
    578 Vertical zooms are defined through the attributs zoom\_begin and zoom\_end of the tag family axis. 
     586Vertical zooms are defined through the attributs zoom\_begin and zoom\_n of the tag family axis. 
    579587It must therefore be done in the axis part of the XML file. 
    580 For example, in \path{NEMOGCM/CONFIG/ORCA2_LIM/iodef_demo.xml}, we provide the following example: 
    581  
    582 \begin{xmllines} 
    583 <axis_group id="deptht" long_name="Vertical T levels" unit="m" positive="down" > 
    584    <axis id="deptht" /> 
    585    <axis id="deptht_myzoom" zoom_begin="1" zoom_end="10" /> 
     588For example, in \path{cfgs/ORCA2_ICE_PISCES/EXPREF/iodef_demo.xml}, we provide the following example: 
     589 
     590\begin{xmllines} 
     591<axis_definition> 
     592   <axis id="deptht" long_name="Vertical T levels" unit="m" positive="down" /> 
     593   <axis id="deptht_zoom" azix_ref="deptht" > 
     594      <zoom_axis zoom_begin="1" zoom_n="10" /> 
     595   </axis> 
    586596\end{xmllines} 
    587597 
    … …  
    595605\end{xmllines} 
    596606 
     607%% ================================================================================================= 
    597608\subsubsection{Control of the output file names} 
    598609 
    … …  
    601612 
    602613\begin{xmllines} 
    603 <file_group id="1d" output_freq="1d" name="myfile_1d" >  
     614<file_group id="1d" output_freq="1d" name="myfile_1d" > 
    604615   <file id="myfileA" name_suffix="_AAA" > <!-- will create file "myfile_1d_AAA"  --> 
    605616      ... 
    … …  
    620631 
    621632\begin{table} 
    622   \scriptsize 
    623633  \begin{tabularx}{\textwidth}{|lX|} 
    624634    \hline 
    … …  
    656666\end{table} 
    657667 
    658 \noindent For example,  
     668\noindent For example, 
    659669\xmlline|<file id="myfile_hzoom" name="myfile_@expname@_@startdate@_freq@freq@" output_freq="1d" >| 
    660670 
    … …  
    668678\noindent will give the following file name radical: \ifile{myfile\_ORCA2\_19891231\_freq1d} 
    669679 
    670 \subsubsection{Other controls of the XML attributes from NEMO} 
    671  
    672 The values of some attributes are defined by subroutine calls within NEMO  
     680%% ================================================================================================= 
     681\subsubsection{Other controls of the XML attributes from \NEMO} 
     682 
     683The values of some attributes are defined by subroutine calls within \NEMO 
    673684(calls to iom\_set\_domain\_attr, iom\_set\_axis\_attr and iom\_set\_field\_attr in \mdl{iom}). 
    674685Any definition given in the XML file will be overwritten. 
    … …  
    681692 
    682693\begin{table} 
    683   \scriptsize 
    684   \begin{tabularx}{\textwidth}{|X|c|c|c|} 
     694  \begin{tabular}{|l|c|c|} 
    685695    \hline 
    686696    tag ids affected by automatic definition of some of their attributes & 
    687697    name attribute                                                       & 
    688     attribute value                      \\ 
     698    attribute value                                                      \\ 
    689699    \hline 
    690700    \hline 
    691701    field\_definition                                                    & 
    692702    freq\_op                                                             & 
    693     \np{rn\_rdt}                         \\ 
     703    \np{rn_rdt}{rn\_rdt}                                                         \\ 
    694704    \hline 
    695705    SBC                                                                  & 
    696706    freq\_op                                                             & 
    697     \np{rn\_rdt} $\times$ \np{nn\_fsbc}  \\ 
     707    \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_fsbc}{nn\_fsbc}                                  \\ 
    698708    \hline 
    699709    ptrc\_T                                                              & 
    700710    freq\_op                                                             & 
    701     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ 
     711    \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_dttrc}{nn\_dttrc}                                \\ 
    702712    \hline 
    703713    diad\_T                                                              & 
    704714    freq\_op                                                             & 
    705     \np{rn\_rdt} $\times$ \np{nn\_dttrc} \\ 
     715    \np{rn_rdt}{rn\_rdt} $\times$ \np{nn_dttrc}{nn\_dttrc}                                \\ 
    706716    \hline 
    707717    EqT, EqU, EqW                                                        & 
    708718    jbegin, ni,                                                          & 
    709     according to the grid                \\ 
    710     & 
     719    according to the grid                                                \\ 
     720                                                                         & 
    711721    name\_suffix                                                         & 
    712     \\ 
     722                                                                         \\ 
    713723    \hline 
    714724    TAO, RAMA and PIRATA moorings                                        & 
    715725    zoom\_ibegin, zoom\_jbegin,                                          & 
    716     according to the grid                \\ 
    717     & 
     726    according to the grid                                                \\ 
     727                                                                         & 
    718728    name\_suffix                                                         & 
    719     \\ 
    720     \hline 
    721   \end{tabularx} 
     729                                                                         \\ 
     730    \hline 
     731  \end{tabular} 
    722732\end{table} 
    723733 
     734%% ================================================================================================= 
    724735\subsubsection{Advanced use of XIOS functionalities} 
    725736 
     737%% ================================================================================================= 
    726738\subsection{XML reference tables} 
    727 \label{subsec:IOM_xmlref} 
     739\label{subsec:DIA_IOM_xmlref} 
    728740 
    729741\begin{enumerate} 
    730 \item 
    731   Simple computation: directly define the computation when refering to the variable in the file definition. 
     742\item Simple computation: directly define the computation when refering to the variable in the file definition. 
    732743 
    733744\begin{xmllines} 
    … …  
    737748\end{xmllines} 
    738749 
    739 \item 
    740   Simple computation: define a new variable and use it in the file definition. 
     750\item Simple computation: define a new variable and use it in the file definition. 
    741751 
    742752in field\_definition: 
    … …  
    755765sst2 won't be evaluated. 
    756766 
    757 \item 
    758   Change of variable precision: 
     767\item Change of variable precision: 
    759768 
    760769\begin{xmllines} 
    … …  
    765774\end{xmllines} 
    766775 
    767 Note that, then the code is crashing, writting real4 variables forces a numerical convection from  
     776Note that, then the code is crashing, writting real4 variables forces a numerical conversion from 
    768777real8 to real4 which will create an internal error in NetCDF and will avoid the creation of the output files. 
    769778Forcing double precision outputs with prec="8" (for example in the field\_definition) will avoid this problem. 
    770779 
    771 \item 
    772   add user defined attributes: 
    773  
    774 \begin{xmllines} 
    775 <file_group id="1d" output_freq="1d" output_level="10" enabled=".true."> <!-- 1d files -->  
     780\item add user defined attributes: 
     781 
     782\begin{xmllines} 
     783<file_group id="1d" output_freq="1d" output_level="10" enabled=".true."> <!-- 1d files --> 
    776784   <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
    777785      <field field_ref="sst" name="tos" > 
    … …  
    782790      <variable id="my_global_attribute" type="string" > blabla_global </variable> 
    783791   </file> 
    784 </file_group>  
    785 \end{xmllines} 
    786  
    787 \item 
    788   use of the ``@'' function: example 1, weighted temporal average 
     792</file_group> 
     793\end{xmllines} 
     794 
     795\item use of the ``@'' function: example 1, weighted temporal average 
    789796 
    790797 - define a new variable in field\_definition 
    … …  
    797804 
    798805\begin{xmllines} 
    799 <file_group id="5d" output_freq="5d"  output_level="10" enabled=".true." >  <!-- 5d files -->   
     806<file_group id="5d" output_freq="5d"  output_level="10" enabled=".true." >  <!-- 5d files --> 
    800807   <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
    801808      <field field_ref="toce" operation="instant" freq_op="5d" > @toce_e3t / @e3t </field> 
    802809   </file> 
    803 </file_group>  
     810</file_group> 
    804811\end{xmllines} 
    805812 
    … …  
    814821Note that in this case, freq\_op must be equal to the file output\_freq. 
    815822 
    816 \item 
    817   use of the ``@'' function: example 2, monthly SSH standard deviation 
     823\item use of the ``@'' function: example 2, monthly SSH standard deviation 
    818824 
    819825 - define a new variable in field\_definition 
    … …  
    826832 
    827833\begin{xmllines} 
    828 <file_group id="1m" output_freq="1m"  output_level="10" enabled=".true." >  <!-- 1m files -->   
     834<file_group id="1m" output_freq="1m"  output_level="10" enabled=".true." >  <!-- 1m files --> 
    829835   <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
    830       <field field_ref="ssh" name="sshstd" long_name="sea_surface_temperature_standard_deviation"  
     836      <field field_ref="ssh" name="sshstd" long_name="sea_surface_temperature_standard_deviation" 
    831837      operation="instant" freq_op="1m" > 
    832838         sqrt( @ssh2 - @ssh * @ssh ) 
    833839      </field> 
    834840   </file> 
    835 </file_group>  
     841</file_group> 
    836842\end{xmllines} 
    837843 
    … …  
    840846here we use the default, average. 
    841847So, in the above case, @ssh2 will do the monthly mean of ssh*ssh. 
    842 Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'':  
     848Operation="instant" refers to the temporal operation to be performed on the field ''sqrt( @ssh2 - @ssh * @ssh )'': 
    843849here the temporal average is alreday done by the ``@'' function so we just use instant. 
    844850field\_ref="ssh" means that attributes not explicitely defined, are inherited from ssh field. 
    845851Note that in this case, freq\_op must be equal to the file output\_freq. 
    846852 
    847 \item 
    848   use of the ``@'' function: example 3, monthly average of SST diurnal cycle 
     853\item use of the ``@'' function: example 3, monthly average of SST diurnal cycle 
    849854 
    850855 - define 2 new variables in field\_definition 
    … …  
    858863 
    859864\begin{xmllines} 
    860 <file_group id="1m" output_freq="1m"  output_level="10" enabled=".true." >  <!-- 1m files -->   
     865<file_group id="1m" output_freq="1m"  output_level="10" enabled=".true." >  <!-- 1m files --> 
    861866   <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
    862867      <field field_ref="sst" name="sstdcy" long_name="amplitude of sst diurnal cycle" operation="average" freq_op="1d" > 
    … …  
    864869      </field> 
    865870   </file> 
    866 </file_group>  
     871</file_group> 
    867872\end{xmllines} 
    868873 
    … …  
    877882field\_ref="sst" means that attributes not explicitely defined, are inherited from sst field. 
    878883 
     884%% ================================================================================================= 
    879885\subsubsection{Tag list per family} 
    880886 
    881887\begin{table} 
    882   \scriptsize 
    883888  \begin{tabularx}{\textwidth}{|l|X|X|l|X|} 
    884889    \hline 
    … …  
    903908    \hline 
    904909  \end{tabularx} 
    905   \caption{Context tags} 
     910  \caption{XIOS: context tags} 
    906911\end{table} 
    907912 
    908913\begin{table} 
    909   \scriptsize 
    910914  \begin{tabularx}{\textwidth}{|l|X|X|X|l|} 
    911915    \hline 
    … …  
    938942    \hline 
    939943  \end{tabularx} 
    940   \caption{Field tags ("\tt{field\_*}")} 
     944  \caption{XIOS: field tags ("\texttt{field\_*}")} 
    941945\end{table} 
    942946 
    943947\begin{table} 
    944   \scriptsize 
    945948  \begin{tabularx}{\textwidth}{|l|X|X|X|l|} 
    946949    \hline 
    … …  
    974977    \hline 
    975978  \end{tabularx} 
    976   \caption{File tags ("\tt{file\_*}")} 
     979  \caption{XIOS: file tags ("\texttt{file\_*}")} 
    977980\end{table} 
    978981 
    979982\begin{table} 
    980   \scriptsize 
    981983  \begin{tabularx}{\textwidth}{|l|X|X|X|X|} 
    982984    \hline 
    … …  
    10071009    \hline 
    10081010  \end{tabularx} 
    1009   \caption{Axis tags ("\tt{axis\_*}")} 
     1011  \caption{XIOS: axis tags ("\texttt{axis\_*}")} 
    10101012\end{table} 
    10111013 
    10121014\begin{table} 
    1013   \scriptsize 
    10141015  \begin{tabularx}{\textwidth}{|l|X|X|X|X|} 
    10151016    \hline 
    … …  
    10401041    \hline 
    10411042  \end{tabularx} 
    1042   \caption{Domain tags ("\tt{domain\_*)}"} 
     1043  \caption{XIOS: domain tags ("\texttt{domain\_*)}"} 
    10431044\end{table} 
    10441045 
    10451046\begin{table} 
    1046   \scriptsize 
    10471047  \begin{tabularx}{\textwidth}{|l|X|X|X|X|} 
    10481048    \hline 
    … …  
    10731073    \hline 
    10741074  \end{tabularx} 
    1075   \caption{Grid tags ("\tt{grid\_*}")} 
     1075  \caption{XIOS: grid tags ("\texttt{grid\_*}")} 
    10761076\end{table} 
    10771077 
     1078%% ================================================================================================= 
    10781079\subsubsection{Attributes list per family} 
    10791080 
    10801081\begin{table} 
    1081   \scriptsize 
    10821082  \begin{tabularx}{\textwidth}{|l|X|l|l|} 
    10831083    \hline 
    … …  
    11141114    \hline 
    11151115  \end{tabularx} 
    1116   \caption{Reference attributes ("\tt{*\_ref}")} 
     1116  \caption{XIOS: reference attributes ("\texttt{*\_ref}")} 
    11171117\end{table} 
    11181118 
    11191119\begin{table} 
    1120   \scriptsize 
    11211120  \begin{tabularx}{\textwidth}{|l|X|l|l|} 
    11221121    \hline 
    … …  
    11501149    \hline 
    11511150  \end{tabularx} 
    1152   \caption{Domain attributes ("\tt{zoom\_*}")} 
     1151  \caption{XIOS: domain attributes ("\texttt{zoom\_*}")} 
    11531152\end{table} 
    11541153 
    11551154\begin{table} 
    1156   \scriptsize 
    11571155  \begin{tabularx}{\textwidth}{|l|X|l|l|} 
    11581156    \hline 
    … …  
    12051203    \hline 
    12061204  \end{tabularx} 
    1207   \caption{File attributes} 
     1205  \caption{XIOS: file attributes} 
    12081206\end{table} 
    12091207 
    12101208\begin{table} 
    1211   \scriptsize 
    12121209  \begin{tabularx}{\textwidth}{|l|X|l|l|} 
    12131210    \hline 
    … …  
    12541251    \hline 
    12551252  \end{tabularx} 
    1256   \caption{Field attributes} 
     1253  \caption{XIOS: field attributes} 
    12571254\end{table} 
    12581255 
    12591256\begin{table} 
    1260   \scriptsize 
    12611257  \begin{tabularx}{\textwidth}{|l|X|X|X|} 
    12621258    \hline 
    … …  
    13131309    \hline 
    13141310  \end{tabularx} 
    1315   \caption{Miscellaneous attributes} 
     1311  \caption{XIOS: miscellaneous attributes} 
    13161312\end{table} 
    13171313 
     1314%% ================================================================================================= 
    13181315\subsection{CF metadata standard compliance} 
    13191316 
    1320 Output from the XIOS-1.0 IO server is compliant with  
     1317Output from the XIOS IO server is compliant with 
    13211318\href{http://cfconventions.org/Data/cf-conventions/cf-conventions-1.5/build/cf-conventions.html}{version 1.5} of 
    1322 the CF metadata standard.  
     1319the CF metadata standard. 
    13231320Therefore while a user may wish to add their own metadata to the output files (as demonstrated in example 4 of 
    1324 section \autoref{subsec:IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
     1321section \autoref{subsec:DIA_IOM_xmlref}) the metadata should, for the most part, comply with the CF-1.5 standard. 
    13251322 
    13261323Some metadata that may significantly increase the file size (horizontal cell areas and vertices) are controlled by 
    1327 the namelist parameter \np{ln\_cfmeta} in the \ngn{namrun} namelist. 
     1324the namelist parameter \np{ln_cfmeta}{ln\_cfmeta} in the \nam{run}{run} namelist. 
    13281325This must be set to true if these metadata are to be included in the output files. 
    13291326 
    1330  
    1331 % ================================================================ 
    1332 %       NetCDF4 support 
    1333 % ================================================================ 
    1334 \section{NetCDF4 support (\protect\key{netcdf4})} 
     1327%% ================================================================================================= 
     1328\section[NetCDF4 support (\texttt{\textbf{key\_netcdf4}})]{NetCDF4 support (\protect\key{netcdf4})} 
    13351329\label{sec:DIA_nc4} 
    13361330 
    … …  
    13401334Chunking and compression can lead to significant reductions in file sizes for a small runtime overhead. 
    13411335For a fuller discussion on chunking and other performance issues the reader is referred to 
    1342 the NetCDF4 documentation found \href{http://www.unidata.ucar.edu/software/netcdf/docs/netcdf.html#Chunking}{here}. 
     1336the NetCDF4 documentation found \href{https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_perf_chunking.html}{here}. 
    13431337 
    13441338The new features are only available when the code has been linked with a NetCDF4 library 
    … …  
    13461340Datasets created with chunking and compression are not backwards compatible with NetCDF3 "classic" format but 
    13471341most analysis codes can be relinked simply with the new libraries and will then read both NetCDF3 and NetCDF4 files. 
    1348 NEMO executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by 
    1349 setting the \np{ln\_nc4zip} logical to false in the \textit{namnc4} namelist: 
    1350  
    1351 %------------------------------------------namnc4---------------------------------------------------- 
    1352  
    1353 \nlst{namnc4} 
    1354 %------------------------------------------------------------------------------------------------------------- 
     1342\NEMO\ executables linked with NetCDF4 libraries can be made to produce NetCDF3 files by 
     1343setting the \np{ln_nc4zip}{ln\_nc4zip} logical to false in the \nam{nc4}{nc4} namelist: 
     1344 
     1345\begin{listing} 
     1346  \nlst{namnc4} 
     1347  \caption{\forcode{&namnc4}} 
     1348  \label{lst:namnc4} 
     1349\end{listing} 
    13551350 
    13561351If \key{netcdf4} has not been defined, these namelist parameters are not read. 
    1357 In this case, \np{ln\_nc4zip} is set false and dummy routines for a few NetCDF4-specific functions are defined. 
     1352In this case, \np{ln_nc4zip}{ln\_nc4zip} is set false and dummy routines for a few NetCDF4-specific functions are defined. 
    13581353These functions will not be used but need to be included so that compilation is possible with NetCDF3 libraries. 
    13591354 
    … …  
    13891384\end{forlines} 
    13901385 
    1391 \noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\tt 46x38x1} respectively in 
    1392 the mono-processor case (\ie global domain of {\small\tt 182x149x31}). 
    1393 An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in  
    1394 table \autoref{tab:NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with 
     1386\noindent for a standard ORCA2\_LIM configuration gives chunksizes of {\small\texttt 46x38x1} respectively in 
     1387the mono-processor case (\ie\ global domain of {\small\texttt 182x149x31}). 
     1388An illustration of the potential space savings that NetCDF4 chunking and compression provides is given in 
     1389table \autoref{tab:DIA_NC4} which compares the results of two short runs of the ORCA2\_LIM reference configuration with 
    13951390a 4x2 mpi partitioning. 
    13961391Note the variation in the compression ratio achieved which reflects chiefly the dry to wet volume ratio of 
    13971392each processing region. 
    13981393 
    1399 %------------------------------------------TABLE---------------------------------------------------- 
    14001394\begin{table} 
    1401   \scriptsize 
    14021395  \centering 
    14031396  \begin{tabular}{lrrr} 
    … …  
    14311424    ORCA2\_2d\_grid\_W\_0007.nc & 4416    & 1368     & 70\%      \\ 
    14321425  \end{tabular} 
    1433   \caption{ 
    1434     \protect\label{tab:NC4} 
    1435     Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression 
    1436   } 
     1426  \caption{Filesize comparison between NetCDF3 and NetCDF4 with chunking and compression} 
     1427  \label{tab:DIA_NC4} 
    14371428\end{table} 
    1438 %---------------------------------------------------------------------------------------------------- 
    14391429 
    14401430When \key{iomput} is activated with \key{netcdf4} chunking and compression parameters for fields produced via 
    1441 \np{iom\_put} calls are set via an equivalent and identically named namelist to \textit{namnc4} in 
    1442 \np{xmlio\_server.def}. 
    1443 Typically this namelist serves the mean files whilst the \ngn{ namnc4} in the main namelist file continues to 
     1431\rou{iom\_put} calls are set via an equivalent and identically named namelist to \nam{nc4}{nc4} in 
     1432\textit{xmlio\_server.def}. 
     1433Typically this namelist serves the mean files whilst the \nam{nc4}{nc4} in the main namelist file continues to 
    14441434serve the restart files. 
    14451435This duplication is unfortunate but appropriate since, if using io\_servers, the domain sizes of 
    14461436the individual files produced by the io\_server processes may be different to those produced by 
    14471437the invidual processing regions and different chunking choices may be desired. 
    1448   
    1449 % ------------------------------------------------------------------------------------------------------------- 
    1450 %       Tracer/Dynamics Trends 
    1451 % ------------------------------------------------------------------------------------------------------------- 
    1452 \section{Tracer/Dynamics trends  (\protect\ngn{namtrd})} 
     1438 
     1439%% ================================================================================================= 
     1440\section[Tracer/Dynamics trends (\forcode{&namtrd})]{Tracer/Dynamics trends (\protect\nam{trd}{trd})} 
    14531441\label{sec:DIA_trd} 
    14541442 
    1455 %------------------------------------------namtrd---------------------------------------------------- 
    1456  
    1457 \nlst{namtrd}  
    1458 %------------------------------------------------------------------------------------------------------------- 
     1443\begin{listing} 
     1444  \nlst{namtrd} 
     1445  \caption{\forcode{&namtrd}} 
     1446  \label{lst:namtrd} 
     1447\end{listing} 
    14591448 
    14601449Each trend of the dynamics and/or temperature and salinity time evolution equations can be send to 
    14611450\mdl{trddyn} and/or \mdl{trdtra} modules (see TRD directory) just after their computation 
    1462 (\ie at the end of each $dyn\cdots.F90$ and/or $tra\cdots.F90$ routines). 
    1463 This capability is controlled by options offered in \ngn{namtrd} namelist. 
    1464 Note that the output are done with xIOS, and therefore the \key{IOM} is required. 
    1465  
    1466 What is done depends on the \ngn{namtrd} logical set to \forcode{.true.}: 
     1451(\ie\ at the end of each \textit{dyn....F90} and/or \textit{tra....F90} routines). 
     1452This capability is controlled by options offered in \nam{trd}{trd} namelist. 
     1453Note that the output are done with XIOS, and therefore the \key{iomput} is required. 
     1454 
     1455What is done depends on the \nam{trd}{trd} logical set to \forcode{.true.}: 
    14671456 
    14681457\begin{description} 
    1469 \item[\np{ln\_glo\_trd}]: 
    1470   at each \np{nn\_trd} time-step a check of the basin averaged properties of 
     1458\item [{\np{ln_glo_trd}{ln\_glo\_trd}}]: at each \np{nn_trd}{nn\_trd} time-step a check of the basin averaged properties of 
    14711459  the momentum and tracer equations is performed. 
    14721460  This also includes a check of $T^2$, $S^2$, $\tfrac{1}{2} (u^2+v2)$, 
    14731461  and potential energy time evolution equations properties; 
    1474 \item[\np{ln\_dyn\_trd}]: 
    1475   each 3D trend of the evolution of the two momentum components is output; 
    1476 \item[\np{ln\_dyn\_mxl}]: 
    1477   each 3D trend of the evolution of the two momentum components averaged over the mixed layer is output; 
    1478 \item[\np{ln\_vor\_trd}]: 
    1479   a vertical summation of the moment tendencies is performed, 
     1462\item [{\np{ln_dyn_trd}{ln\_dyn\_trd}}]: each 3D trend of the evolution of the two momentum components is output; 
     1463\item [{\np{ln_dyn_mxl}{ln\_dyn\_mxl}}]: each 3D trend of the evolution of the two momentum components averaged over the mixed layer is output; 
     1464\item [{\np{ln_vor_trd}{ln\_vor\_trd}}]: a vertical summation of the moment tendencies is performed, 
    14801465  then the curl is computed to obtain the barotropic vorticity tendencies which are output; 
    1481 \item[\np{ln\_KE\_trd}] : 
    1482   each 3D trend of the Kinetic Energy equation is output; 
    1483 \item[\np{ln\_tra\_trd}]: 
    1484   each 3D trend of the evolution of temperature and salinity is output; 
    1485 \item[\np{ln\_tra\_mxl}]: 
    1486   each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output; 
     1466\item [{\np{ln_KE_trd}{ln\_KE\_trd}}]  : each 3D trend of the Kinetic Energy equation is output; 
     1467\item [{\np{ln_tra_trd}{ln\_tra\_trd}}]: each 3D trend of the evolution of temperature and salinity is output; 
     1468\item [{\np{ln_tra_mxl}{ln\_tra\_mxl}}]: each 2D trend of the evolution of temperature and salinity averaged over the mixed layer is output; 
    14871469\end{description} 
    14881470 
    1489 Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via  
    1490 the \key{trdtrc} and \key{trdmld\_trc} CPP keys. 
     1471Note that the mixed layer tendency diagnostic can also be used on biogeochemical models via 
     1472the \key{trdtrc} and \key{trdmxl\_trc} CPP keys. 
    14911473 
    14921474\textbf{Note that} in the current version (v3.6), many changes has been introduced but not fully tested. 
    1493 In particular, options associated with \np{ln\_dyn\_mxl}, \np{ln\_vor\_trd}, and \np{ln\_tra\_mxl} are not working, 
    1494 and none of the options have been tested with variable volume (\ie \key{vvl} defined). 
    1495  
    1496 % ------------------------------------------------------------------------------------------------------------- 
    1497 %       On-line Floats trajectories 
    1498 % ------------------------------------------------------------------------------------------------------------- 
    1499 \section{FLO: On-Line Floats trajectories (\protect\key{floats})} 
    1500 \label{sec:FLO} 
    1501 %--------------------------------------------namflo------------------------------------------------------- 
    1502  
    1503 \nlst{namflo}  
    1504 %-------------------------------------------------------------------------------------------------------------- 
     1475In particular, options associated with \np{ln_dyn_mxl}{ln\_dyn\_mxl}, \np{ln_vor_trd}{ln\_vor\_trd}, and \np{ln_tra_mxl}{ln\_tra\_mxl} are not working, 
     1476and none of the options have been tested with variable volume (\ie\ \np[=.true.]{ln_linssh}{ln\_linssh}). 
     1477 
     1478%% ================================================================================================= 
     1479\section[FLO: On-Line Floats trajectories (\texttt{\textbf{key\_floats}})]{FLO: On-Line Floats trajectories (\protect\key{floats})} 
     1480\label{sec:DIA_FLO} 
     1481 
     1482\begin{listing} 
     1483  \nlst{namflo} 
     1484  \caption{\forcode{&namflo}} 
     1485  \label{lst:namflo} 
     1486\end{listing} 
    15051487 
    15061488The on-line computation of floats advected either by the three dimensional velocity field or constraint to 
    15071489remain at a given depth ($w = 0$ in the computation) have been introduced in the system during the CLIPPER project. 
    1508 Options are defined by \ngn{namflo} namelis variables. 
    1509 The algorithm used is based either on the work of \cite{Blanke_Raynaud_JPO97} (default option), 
    1510 or on a $4^th$ Runge-Hutta algorithm (\np{ln\_flork4}\forcode{ = .true.}). 
    1511 Note that the \cite{Blanke_Raynaud_JPO97} algorithm have the advantage of providing trajectories which 
     1490Options are defined by \nam{flo}{flo} namelist variables. 
     1491The algorithm used is based either on the work of \cite{blanke.raynaud_JPO97} (default option), 
     1492or on a $4^th$ Runge-Hutta algorithm (\np[=.true.]{ln_flork4}{ln\_flork4}). 
     1493Note that the \cite{blanke.raynaud_JPO97} algorithm have the advantage of providing trajectories which 
    15121494are consistent with the numeric of the code, so that the trajectories never intercept the bathymetry. 
    15131495 
     1496%% ================================================================================================= 
    15141497\subsubsection{Input data: initial coordinates} 
    15151498 
    15161499Initial coordinates can be given with Ariane Tools convention 
    1517 (IJK coordinates, (\np{ln\_ariane}\forcode{ = .true.}) ) or with longitude and latitude. 
    1518  
    1519 In case of Ariane convention, input filename is \np{init\_float\_ariane}. 
     1500(IJK coordinates, (\np[=.true.]{ln_ariane}{ln\_ariane}) ) or with longitude and latitude. 
     1501 
     1502In case of Ariane convention, input filename is \textit{init\_float\_ariane}. 
    15201503Its format is: \\ 
    1521 {\scriptsize \texttt{I J K nisobfl itrash itrash}} 
     1504{ \texttt{I J K nisobfl itrash}} 
    15221505 
    15231506\noindent with: 
    … …  
    15251508 - I,J,K  : indexes of initial position 
    15261509 
    1527  - nisobfl: 0 for an isobar float, 1 for a float following the w velocity   
     1510 - nisobfl: 0 for an isobar float, 1 for a float following the w velocity 
    15281511 
    15291512 - itrash : set to zero; it is a dummy variable to respect Ariane Tools convention 
    … …  
    15311514\noindent Example: \\ 
    15321515\noindent 
    1533 {\scriptsize 
     1516{ 
    15341517  \texttt{ 
    15351518    100.00000  90.00000  -1.50000 1.00000   0.00000   \\ 
    … …  
    15421525In the other case (longitude and latitude), input filename is init\_float. 
    15431526Its format is: \\ 
    1544 {\scriptsize \texttt{Long Lat depth nisobfl ngrpfl itrash}} 
     1527{ \texttt{Long Lat depth nisobfl ngrpfl itrash}} 
    15451528 
    15461529\noindent with: 
    … …  
    15561539\noindent Example: \\ 
    15571540\noindent 
    1558 {\scriptsize 
     1541{ 
    15591542  \texttt{ 
    15601543    20.0 0.0 0.0 0 1 1    \\ 
    … …  
    15651548} \\ 
    15661549 
    1567 \np{jpnfl} is the total number of floats during the run. 
    1568 When initial positions are read in a restart file (\np{ln\_rstflo}\forcode{ = .true.} ), 
    1569 \np{jpnflnewflo} can be added in the initialization file. 
    1570  
     1550\np{jpnfl}{jpnfl} is the total number of floats during the run. 
     1551When initial positions are read in a restart file (\np[=.true.]{ln_rstflo}{ln\_rstflo} ), 
     1552\np{jpnflnewflo}{jpnflnewflo} can be added in the initialization file. 
     1553 
     1554%% ================================================================================================= 
    15711555\subsubsection{Output data} 
    15721556 
    1573 \np{nn\_writefl} is the frequency of writing in float output file and \np{nn\_stockfl} is the frequency of 
     1557\np{nn_writefl}{nn\_writefl} is the frequency of writing in float output file and \np{nn_stockfl}{nn\_stockfl} is the frequency of 
    15741558creation of the float restart file. 
    15751559 
    1576 Output data can be written in ascii files (\np{ln\_flo\_ascii}\forcode{ = .true.}). 
     1560Output data can be written in ascii files (\np[=.true.]{ln_flo_ascii}{ln\_flo\_ascii}). 
    15771561In that case, output filename is trajec\_float. 
    15781562 
    1579 Another possiblity of writing format is Netcdf (\np{ln\_flo\_ascii}\forcode{ = .false.}). 
    1580 There are 2 possibilities: 
    1581  
    1582 - if (\key{iomput}) is used, outputs are selected in  iodef.xml. 
     1563Another possiblity of writing format is Netcdf (\np[=.false.]{ln_flo_ascii}{ln\_flo\_ascii}) with 
     1564\key{iomput} and outputs selected in iodef.xml. 
    15831565Here it is an example of specification to put in files description section: 
    15841566 
    … …  
    15971579\end{xmllines} 
    15981580 
    1599  -  if (\key{iomput}) is not used, a file called \ifile{trajec\_float} will be created by IOIPSL library. 
    1600  
    1601  See also \href{http://stockage.univ-brest.fr/~grima/Ariane/}{here} the web site describing the off-line use of 
    1602  this marvellous diagnostic tool. 
    1603  
    1604 % ------------------------------------------------------------------------------------------------------------- 
    1605 %       Harmonic analysis of tidal constituents 
    1606 % ------------------------------------------------------------------------------------------------------------- 
    1607 \section{Harmonic analysis of tidal constituents (\protect\key{diaharm}) } 
     1581%% ================================================================================================= 
     1582\section[Harmonic analysis of tidal constituents (\texttt{\textbf{key\_diaharm}})]{Harmonic analysis of tidal constituents (\protect\key{diaharm})} 
    16081583\label{sec:DIA_diag_harm} 
    16091584 
    1610 %------------------------------------------namdia_harm---------------------------------------------------- 
    1611 % 
    1612 \nlst{nam_diaharm} 
    1613 %---------------------------------------------------------------------------------------------------------- 
     1585\begin{listing} 
     1586  \nlst{nam_diaharm} 
     1587  \caption{\forcode{&nam_diaharm}} 
     1588  \label{lst:nam_diaharm} 
     1589\end{listing} 
    16141590 
    16151591A module is available to compute the amplitude and phase of tidal waves. 
    16161592This on-line Harmonic analysis is actived with \key{diaharm}. 
    16171593 
    1618 Some parameters are available in namelist \ngn{namdia\_harm}: 
    1619  
    1620  - \np{nit000\_han} is the first time step used for harmonic analysis 
    1621  
    1622  - \np{nitend\_han} is the  last time step used for harmonic analysis 
    1623  
    1624  - \np{nstep\_han}  is the  time step frequency for harmonic analysis 
    1625  
    1626  - \np{nb\_ana}     is the number of harmonics to analyse 
    1627  
    1628  - \np{tname}       is an array with names of tidal constituents to analyse 
    1629  
    1630  \np{nit000\_han} and \np{nitend\_han} must be between \np{nit000} and \np{nitend} of the simulation. 
     1594Some parameters are available in namelist \nam{_diaharm}{\_diaharm}: 
     1595 
     1596 - \np{nit000_han}{nit000\_han} is the first time step used for harmonic analysis 
     1597 
     1598 - \np{nitend_han}{nitend\_han} is the  last time step used for harmonic analysis 
     1599 
     1600 - \np{nstep_han}{nstep\_han}  is the  time step frequency for harmonic analysis 
     1601 
     1602% - \np{nb_ana}{nb\_ana}     is the number of harmonics to analyse 
     1603 
     1604 - \np{tname}{tname}       is an array with names of tidal constituents to analyse 
     1605 
     1606 \np{nit000_han}{nit000\_han} and \np{nitend_han}{nitend\_han} must be between \np{nit000}{nit000} and \np{nitend}{nitend} of the simulation. 
    16311607 The restart capability is not implemented. 
    16321608 
    … …  
    16491625We obtain in output $C_{j}$ and $S_{j}$ for each tidal wave. 
    16501626 
    1651 % ------------------------------------------------------------------------------------------------------------- 
    1652 %       Sections transports 
    1653 % ------------------------------------------------------------------------------------------------------------- 
    1654 \section{Transports across sections (\protect\key{diadct}) } 
     1627%% ================================================================================================= 
     1628\section[Transports across sections (\texttt{\textbf{key\_diadct}})]{Transports across sections (\protect\key{diadct})} 
    16551629\label{sec:DIA_diag_dct} 
    16561630 
    1657 %------------------------------------------namdct---------------------------------------------------- 
    1658  
    1659 \nlst{namdct} 
    1660 %------------------------------------------------------------------------------------------------------------- 
     1631\begin{listing} 
     1632  \nlst{nam_diadct} 
     1633  \caption{\forcode{&nam_diadct}} 
     1634  \label{lst:nam_diadct} 
     1635\end{listing} 
    16611636 
    16621637A module is available to compute the transport of volume, heat and salt through sections. 
    … …  
    16641639 
    16651640Each section is defined by the coordinates of its 2 extremities. 
    1666 The pathways between them are contructed using tools which can be found in \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT} 
    1667 and are written in a binary file \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is later read in by 
    1668 NEMO to compute on-line transports. 
     1641The pathways between them are contructed using tools which can be found in \texttt{tools/SECTIONS\_DIADCT} 
     1642and are written in a binary file \texttt{section\_ijglobal.diadct} which is later read in by 
     1643\NEMO\ to compute on-line transports. 
    16691644 
    16701645The on-line transports module creates three output ascii files: 
    … …  
    16761651- \texttt{salt\_transport}   for   salt transports (unit: $10^{9}Kg s^{-1}$) \\ 
    16771652 
    1678 Namelist variables in \ngn{namdct} control how frequently the flows are summed and the time scales over which 
     1653Namelist variables in \nam{_diadct}{\_diadct} control how frequently the flows are summed and the time scales over which 
    16791654they are averaged, as well as the level of output for debugging: 
    1680 \np{nn\_dct}   : frequency of instantaneous transports computing 
    1681 \np{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) 
    1682 \np{nn\_debug} : debugging of the section 
    1683  
     1655\np{nn_dct}{nn\_dct}   : frequency of instantaneous transports computing 
     1656\np{nn_dctwri}{nn\_dctwri}: frequency of writing ( mean of instantaneous transports ) 
     1657\np{nn_debug}{nn\_debug} : debugging of the section 
     1658 
     1659%% ================================================================================================= 
    16841660\subsubsection{Creating a binary file containing the pathway of each section} 
    16851661 
    1686 In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, 
     1662In \texttt{tools/SECTIONS\_DIADCT/run}, 
    16871663the file \textit{ {list\_sections.ascii\_global}} contains a list of all the sections that are to be computed 
    16881664(this list of sections is based on MERSEA project metrics). 
    … …  
    16911667 
    16921668Each section is defined by: \\ 
    1693 \noindent {\scriptsize \texttt{long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name}} \\ 
     1669\noindent { \texttt{long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name}} \\ 
    16941670with: 
    16951671 
    … …  
    16981674 - \texttt{long2 lat2}, coordinates of the second extremity of the section; 
    16991675 
    1700  - \texttt{nclass}    the number of bounds of your classes (\eg bounds for 2 classes); 
     1676 - \texttt{nclass}    the number of bounds of your classes (\eg\ bounds for 2 classes); 
    17011677 
    17021678 - \texttt{okstrpond} to compute    heat and       salt transports, \texttt{nostrpond} if no; 
    … …  
    17081684 
    17091685\noindent If nclass $\neq$ 0, the next lines contain the class type and the nclass bounds: \\ 
    1710 {\scriptsize 
     1686{ 
    17111687  \texttt{ 
    17121688    long1 lat1 long2 lat2 nclass (ok/no)strpond (no)ice section\_name \\ 
    … …  
    17311707 
    17321708 - \texttt{zsigp} for potential density classes \\ 
    1733    
     1709 
    17341710 The script \texttt{job.ksh} computes the pathway for each section and creates a binary file 
    1735  \texttt{section\_ijglobal.diadct\_ORCA2\_LIM} which is read by NEMO. \\ 
     1711 \texttt{section\_ijglobal.diadct} which is read by \NEMO. \\ 
    17361712 
    17371713 It is possible to use this tools for new configuations: \texttt{job.ksh} has to be updated with 
    … …  
    17411717 and the ATL\_Cuba\_Florida with 4 temperature clases (5 class bounds), are shown: \\ 
    17421718 \noindent 
    1743  {\scriptsize 
     1719 { 
    17441720   \texttt{ 
    17451721     -68.    -54.5   -60.    -64.7  00 okstrpond noice ACC\_Drake\_Passage \\ 
    … …  
    17531729 } 
    17541730 
     1731%% ================================================================================================= 
    17551732\subsubsection{To read the output files} 
    17561733 
    17571734The output format is: \\ 
    1758 {\scriptsize 
     1735{ 
    17591736  \texttt{ 
    17601737    date, time-step number, section number,                \\ 
    … …  
    17781755 
    17791756\begin{table} 
    1780   \scriptsize 
    17811757  \begin{tabular}{|l|l|l|l|l|} 
    17821758    \hline 
    … …  
    17921768\end{table} 
    17931769 
    1794 % ================================================================ 
    1795 % Steric effect in sea surface height 
    1796 % ================================================================ 
     1770%% ================================================================================================= 
    17971771\section{Diagnosing the steric effect in sea surface height} 
    17981772\label{sec:DIA_steric} 
    1799  
    18001773 
    18011774Changes in steric sea level are caused when changes in the density of the water column imply an expansion or 
    … …  
    18091782The steric effect is therefore not explicitely represented. 
    18101783This approximation does not represent a serious error with respect to the flow field calculated by the model 
    1811 \citep{Greatbatch_JGR94}, but extra attention is required when investigating sea level, 
     1784\citep{greatbatch_JGR94}, but extra attention is required when investigating sea level, 
    18121785as steric changes are an important contribution to local changes in sea level on seasonal and climatic time scales. 
    18131786This is especially true for investigation into sea level rise due to global warming. 
    18141787 
    18151788Fortunately, the steric contribution to the sea level consists of a spatially uniform component that 
    1816 can be diagnosed by considering the mass budget of the world ocean \citep{Greatbatch_JGR94}. 
     1789can be diagnosed by considering the mass budget of the world ocean \citep{greatbatch_JGR94}. 
    18171790In order to better understand how global mean sea level evolves and thus how the steric sea level can be diagnosed, 
    18181791we compare, in the following, the non-Boussinesq and Boussinesq cases. 
    18191792 
    18201793Let denote 
    1821 $\mathcal{M}$ the total mass    of liquid seawater ($\mathcal{M} = \int_D \rho dv$),  
    1822 $\mathcal{V}$ the total volume  of        seawater      ($\mathcal{V} = \int_D dv$),  
    1823 $\mathcal{A}$ the total surface of       the ocean      ($\mathcal{A} = \int_S ds$),  
    1824 $\bar{\rho}$ the global mean  seawater (\textit{in situ}) density  
     1794$\mathcal{M}$ the total mass    of liquid seawater ($\mathcal{M} = \int_D \rho dv$), 
     1795$\mathcal{V}$ the total volume  of        seawater      ($\mathcal{V} = \int_D dv$), 
     1796$\mathcal{A}$ the total surface of       the ocean      ($\mathcal{A} = \int_S ds$), 
     1797$\bar{\rho}$ the global mean  seawater (\textit{in situ}) density 
    18251798($\bar{\rho} = 1/\mathcal{V} \int_D \rho \,dv$), and 
    1826 $\bar{\eta}$ the global mean sea level  
     1799$\bar{\eta}$ the global mean sea level 
    18271800($\bar{\eta} = 1/\mathcal{A} \int_S \eta \,ds$). 
    18281801 
    … …  
    18341807    \mathcal{V} &=  \mathcal{A}  \;\bar{\eta} 
    18351808  \end{split} 
    1836   \label{eq:MV_nBq} 
     1809  \label{eq:DIA_MV_nBq} 
    18371810\end{equation} 
    18381811 
    … …  
    18421815  \frac{1}{e_3} \partial_t ( e_3\,\rho) + \nabla( \rho \, \textbf{U} ) 
    18431816  = \left. \frac{\textit{emp}}{e_3}\right|_\textit{surface} 
    1844   \label{eq:Co_nBq} 
     1817  \label{eq:DIA_Co_nBq} 
    18451818\end{equation} 
    18461819 
    18471820where $\rho$ is the \textit{in situ} density, and \textit{emp} the surface mass exchanges with the other media of 
    18481821the Earth system (atmosphere, sea-ice, land). 
    1849 Its global averaged leads to the total mass change  
     1822Its global averaged leads to the total mass change 
    18501823 
    18511824\begin{equation} 
    18521825  \partial_t \mathcal{M} = \mathcal{A} \;\overline{\textit{emp}} 
    1853   \label{eq:Mass_nBq} 
     1826  \label{eq:DIA_Mass_nBq} 
    18541827\end{equation} 
    18551828 
    18561829where $\overline{\textit{emp}} = \int_S \textit{emp}\,ds$ is the net mass flux through the ocean surface. 
    1857 Bringing \autoref{eq:Mass_nBq} and the time derivative of \autoref{eq:MV_nBq} together leads to 
     1830Bringing \autoref{eq:DIA_Mass_nBq} and the time derivative of \autoref{eq:DIA_MV_nBq} together leads to 
    18581831the evolution equation of the mean sea level 
    18591832 
    … …  
    18611834  \partial_t \bar{\eta} = \frac{\overline{\textit{emp}}}{ \bar{\rho}} 
    18621835  - \frac{\mathcal{V}}{\mathcal{A}} \;\frac{\partial_t \bar{\rho} }{\bar{\rho}} 
    1863   \label{eq:ssh_nBq} 
     1836  \label{eq:DIA_ssh_nBq} 
    18641837\end{equation} 
    18651838 
    1866 The first term in equation \autoref{eq:ssh_nBq} alters sea level by adding or subtracting mass from the ocean.  
    1867 The second term arises from temporal changes in the global mean density; \ie from steric effects. 
     1839The first term in equation \autoref{eq:DIA_ssh_nBq} alters sea level by adding or subtracting mass from the ocean. 
     1840The second term arises from temporal changes in the global mean density; \ie\ from steric effects. 
    18681841 
    18691842In a Boussinesq fluid, $\rho$ is replaced by $\rho_o$ in all the equation except when $\rho$ appears multiplied by 
    1870 the gravity (\ie in the hydrostatic balance of the primitive Equations). 
    1871 In particular, the mass conservation equation, \autoref{eq:Co_nBq}, degenerates into the incompressibility equation: 
     1843the gravity (\ie\ in the hydrostatic balance of the primitive Equations). 
     1844In particular, the mass conservation equation, \autoref{eq:DIA_Co_nBq}, degenerates into the incompressibility equation: 
    18721845 
    18731846\[ 
    18741847  \frac{1}{e_3} \partial_t ( e_3 ) + \nabla( \textbf{U} ) = \left. \frac{\textit{emp}}{\rho_o \,e_3}\right|_ \textit{surface} 
    1875   % \label{eq:Co_Bq} 
     1848  % \label{eq:DIA_Co_Bq} 
    18761849\] 
    18771850 
    … …  
    18801853\[ 
    18811854  \partial_t \mathcal{V} = \mathcal{A} \;\frac{\overline{\textit{emp}}}{\rho_o} 
    1882   % \label{eq:V_Bq} 
     1855  % \label{eq:DIA_V_Bq} 
    18831856\] 
    18841857 
    … …  
    18881861the ocean surface, not by changes in mean mass of the ocean: the steric effect is missing in a Boussinesq fluid. 
    18891862 
    1890 Nevertheless, following \citep{Greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 
    1891 considering the mass budget of the ocean.  
     1863Nevertheless, following \citep{greatbatch_JGR94}, the steric effect on the volume can be diagnosed by 
     1864considering the mass budget of the ocean. 
    18921865The apparent changes in $\mathcal{M}$, mass of the ocean, which are not induced by surface mass flux 
    18931866must be compensated by a spatially uniform change in the mean sea level due to expansion/contraction of the ocean 
    1894 \citep{Greatbatch_JGR94}. 
     1867\citep{greatbatch_JGR94}. 
    18951868In others words, the Boussinesq mass, $\mathcal{M}_o$, can be related to $\mathcal{M}$, 
    18961869the total mass of the ocean seen by the Boussinesq model, via the steric contribution to the sea level, 
    … …  
    18991872\begin{equation} 
    19001873  \mathcal{M}_o = \mathcal{M} + \rho_o \,\eta_s \,\mathcal{A} 
    1901   \label{eq:M_Bq} 
     1874  \label{eq:DIA_M_Bq} 
    19021875\end{equation} 
    19031876 
    … …  
    19051878is converted into a mean change in sea level. 
    19061879Introducing the total density anomaly, $\mathcal{D}= \int_D d_a \,dv$, 
    1907 where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO (cf. \autoref{subsec:TRA_eos}) 
    1908 in \autoref{eq:M_Bq} leads to a very simple form for the steric height: 
     1880where $d_a = (\rho -\rho_o ) / \rho_o$ is the density anomaly used in \NEMO\ (cf. \autoref{subsec:TRA_eos}) 
     1881in \autoref{eq:DIA_M_Bq} leads to a very simple form for the steric height: 
    19091882 
    19101883\begin{equation} 
    19111884  \eta_s = - \frac{1}{\mathcal{A}} \mathcal{D} 
    1912   \label{eq:steric_Bq} 
     1885  \label{eq:DIA_steric_Bq} 
    19131886\end{equation} 
    19141887 
    19151888The above formulation of the steric height of a Boussinesq ocean requires four remarks. 
    19161889First, one can be tempted to define $\rho_o$ as the initial value of $\mathcal{M}/\mathcal{V}$, 
    1917 \ie set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. 
     1890\ie\ set $\mathcal{D}_{t=0}=0$, so that the initial steric height is zero. 
    19181891We do not recommend that. 
    19191892Indeed, in this case $\rho_o$ depends on the initial state of the ocean. 
    … …  
    19241897This value is a sensible choice for the reference density used in a Boussinesq ocean climate model since, 
    19251898with the exception of only a small percentage of the ocean, density in the World Ocean varies by no more than 
    1926 2$\%$ from this value (\cite{Gill1982}, page 47). 
     18992$\%$ from this value (\cite{gill_bk82}, page 47). 
    19271900 
    19281901Second, we have assumed here that the total ocean surface, $\mathcal{A}$, 
    19291902does not change when the sea level is changing as it is the case in all global ocean GCMs 
    19301903(wetting and drying of grid point is not allowed). 
    1931    
    1932 Third, the discretisation of \autoref{eq:steric_Bq} depends on the type of free surface which is considered. 
    1933 In the non linear free surface case, \ie \key{vvl} defined, it is given by 
     1904 
     1905Third, the discretisation of \autoref{eq:DIA_steric_Bq} depends on the type of free surface which is considered. 
     1906In the non linear free surface case, \ie\ \np[=.true.]{ln_linssh}{ln\_linssh}, it is given by 
    19341907 
    19351908\[ 
    19361909  \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t} e_{2t} e_{3t} }{ \sum_{i,\,j,\,k}       e_{1t} e_{2t} e_{3t} } 
    1937   % \label{eq:discrete_steric_Bq_nfs} 
     1910  % \label{eq:DIA_discrete_steric_Bq_nfs} 
    19381911\] 
    19391912 
    … …  
    19451918  \eta_s = - \frac{ \sum_{i,\,j,\,k} d_a\; e_{1t}e_{2t}e_{3t} + \sum_{i,\,j} d_a\; e_{1t}e_{2t} \eta } 
    19461919                  { \sum_{i,\,j,\,k}       e_{1t}e_{2t}e_{3t} + \sum_{i,\,j}       e_{1t}e_{2t} \eta } 
    1947   % \label{eq:discrete_steric_Bq_fs} 
     1920  % \label{eq:DIA_discrete_steric_Bq_fs} 
    19481921\] 
    19491922 
    … …  
    19541927so that there are no associated ocean currents. 
    19551928Hence, the dynamically relevant sea level is the effective sea level, 
    1956 \ie the sea level as if sea ice (and snow) were converted to liquid seawater \citep{Campin_al_OM08}. 
    1957 However, in the current version of \NEMO the sea-ice is levitating above the ocean without mass exchanges between 
     1929\ie\ the sea level as if sea ice (and snow) were converted to liquid seawater \citep{campin.marshall.ea_OM08}. 
     1930However, in the current version of \NEMO\ the sea-ice is levitating above the ocean without mass exchanges between 
    19581931ice and ocean. 
    19591932Therefore the model effective sea level is always given by $\eta + \eta_s$, whether or not there is sea ice present. 
    … …  
    19651938\[ 
    19661939  \eta_s = - \frac{1}{\mathcal{A}} \int_D d_a(T,S_o,p_o) \,dv 
    1967   % \label{eq:thermosteric_Bq} 
     1940  % \label{eq:DIA_thermosteric_Bq} 
    19681941\] 
    19691942 
    19701943where $S_o$ and $p_o$ are the initial salinity and pressure, respectively. 
    19711944 
    1972 Both steric and thermosteric sea level are computed in \mdl{diaar5} which needs the \key{diaar5} defined to 
    1973 be called. 
    1974  
    1975 % ------------------------------------------------------------------------------------------------------------- 
    1976 %       Other Diagnostics 
    1977 % ------------------------------------------------------------------------------------------------------------- 
    1978 \section{Other diagnostics (\protect\key{diahth}, \protect\key{diaar5})} 
     1945Both steric and thermosteric sea level are computed in \mdl{diaar5}. 
     1946 
     1947%% ================================================================================================= 
     1948\section{Other diagnostics} 
    19791949\label{sec:DIA_diag_others} 
    19801950 
    … …  
    19821952The available ready-to-add diagnostics modules can be found in directory DIA. 
    19831953 
    1984 \subsection{Depth of various quantities (\protect\mdl{diahth})} 
     1954%% ================================================================================================= 
     1955\subsection[Depth of various quantities (\textit{diahth.F90})]{Depth of various quantities (\protect\mdl{diahth})} 
    19851956 
    19861957Among the available diagnostics the following ones are obtained when defining the \key{diahth} CPP key: 
    19871958 
    1988 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
     1959- the mixed layer depth (based on a density criterion \citep{de-boyer-montegut.madec.ea_JGR04}) (\mdl{diahth}) 
    19891960 
    19901961- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) 
    … …  
    19941965- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) 
    19951966 
    1996 % ----------------------------------------------------------- 
    1997 %     Poleward heat and salt transports 
    1998 % ----------------------------------------------------------- 
    1999  
    2000 \subsection{Poleward heat and salt transports (\protect\mdl{diaptr})} 
    2001  
    2002 %------------------------------------------namptr----------------------------------------- 
    2003  
    2004 \nlst{namptr}  
    2005 %----------------------------------------------------------------------------------------- 
    2006  
    2007 The poleward heat and salt transports, their advective and diffusive component, 
    2008 and the meriodional stream function can be computed on-line in \mdl{diaptr} \np{ln\_diaptr} to true 
    2009 (see the \textit{\ngn{namptr} } namelist below). 
    2010 When \np{ln\_subbas}\forcode{ = .true.}, transports and stream function are computed for the Atlantic, Indian, 
     1967\begin{figure}[!t] 
     1968  \centering 
     1969  \includegraphics[width=0.66\textwidth]{DIA_mask_subasins} 
     1970  \caption[Decomposition of the World Ocean to compute transports as well as 
     1971  the meridional stream-function]{ 
     1972    Decomposition of the World Ocean (here ORCA2) into sub-basin used in to 
     1973    compute the heat and salt transports as well as the meridional stream-function: 
     1974    Atlantic basin (red), Pacific basin (green), 
     1975    Indian basin (blue), Indo-Pacific basin (blue+green). 
     1976    Note that semi-enclosed seas (Red, Med and Baltic seas) as well as 
     1977    Hudson Bay are removed from the sub-basins. 
     1978    Note also that the Arctic Ocean has been split into Atlantic and 
     1979    Pacific basins along the North fold line. 
     1980  } 
     1981  \label{fig:DIA_mask_subasins} 
     1982\end{figure} 
     1983 
     1984%% ================================================================================================= 
     1985\subsection[CMIP specific diagnostics (\textit{diaar5.F90}, \textit{diaptr.F90})]{CMIP specific diagnostics (\protect\mdl{diaar5})} 
     1986 
     1987A series of diagnostics has been added in the \mdl{diaar5} and \mdl{diaptr}. 
     1988In \mdl{diaar5} they correspond to outputs that are required for AR5 simulations (CMIP5) 
     1989(see also \autoref{sec:DIA_steric} for one of them). 
     1990The module \mdl{diaar5} is active when one of the following outputs is required : 
     1991global total volume (voltot), global mean ssh (sshtot), global total mass (masstot), global mean temperature (temptot), 
     1992global mean ssh steric (sshsteric), global mean ssh thermosteric (sshthster), global mean salinity (saltot), 
     1993sea water pressure at sea floor (botpres), dynamic sea surface height (sshdyn). 
     1994 
     1995In \mdl{diaptr} when \np[=.true.]{ln_diaptr}{ln\_diaptr} 
     1996(see the \nam{ptr}{ptr} namelist below) can be computed on-line the poleward heat and salt transports, 
     1997their advective and diffusive component, and the meriodional stream function . 
     1998When \np[=.true.]{ln_subbas}{ln\_subbas}, transports and stream function are computed for the Atlantic, Indian, 
    20111999Pacific and Indo-Pacific Oceans (defined north of 30\deg{S}) as well as for the World Ocean. 
    20122000The sub-basin decomposition requires an input file (\ifile{subbasins}) which contains three 2D mask arrays, 
    2013 the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:mask_subasins}). 
    2014  
    2015 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2016 \begin{figure}[!t] 
    2017   \begin{center} 
    2018     \includegraphics[width=1.0\textwidth]{Fig_mask_subasins} 
    2019     \caption{ 
    2020       \protect\label{fig:mask_subasins} 
    2021       Decomposition of the World Ocean (here ORCA2) into sub-basin used in to 
    2022       compute the heat and salt transports as well as the meridional stream-function: 
    2023       Atlantic basin (red), Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
    2024       Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay are removed from the sub-basins. 
    2025       Note also that the Arctic Ocean has been split into Atlantic and Pacific basins along the North fold line. 
    2026     } 
    2027   \end{center} 
    2028 \end{figure}   
    2029 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2030  
    2031 % ----------------------------------------------------------- 
    2032 %       CMIP specific diagnostics  
    2033 % ----------------------------------------------------------- 
    2034 \subsection{CMIP specific diagnostics (\protect\mdl{diaar5})} 
    2035  
    2036 A series of diagnostics has been added in the \mdl{diaar5}. 
    2037 They corresponds to outputs that are required for AR5 simulations (CMIP5) 
    2038 (see also \autoref{sec:DIA_steric} for one of them). 
    2039 Activating those outputs requires to define the \key{diaar5} CPP key. 
    2040  
    2041 % ----------------------------------------------------------- 
    2042 %       25 hour mean and hourly Surface, Mid and Bed  
    2043 % ----------------------------------------------------------- 
     2001the Indo-Pacific mask been deduced from the sum of the Indian and Pacific mask (\autoref{fig:DIA_mask_subasins}). 
     2002 
     2003\begin{listing} 
     2004  \nlst{namptr} 
     2005  \caption{\forcode{&namptr}} 
     2006  \label{lst:namptr} 
     2007\end{listing} 
     2008 
     2009%% ================================================================================================= 
    20442010\subsection{25 hour mean output for tidal models} 
    20452011 
    2046 %------------------------------------------nam_dia25h------------------------------------- 
    2047  
    2048 \nlst{nam_dia25h} 
    2049 %----------------------------------------------------------------------------------------- 
     2012\begin{listing} 
     2013  \nlst{nam_dia25h} 
     2014  \caption{\forcode{&nam_dia25h}} 
     2015  \label{lst:nam_dia25h} 
     2016\end{listing} 
    20502017 
    20512018A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean. 
    … …  
    20542021This diagnostic is actived with the logical $ln\_dia25h$. 
    20552022 
    2056 % ----------------------------------------------------------- 
    2057 %     Top Middle and Bed hourly output 
    2058 % ----------------------------------------------------------- 
     2023%% ================================================================================================= 
    20592024\subsection{Top middle and bed hourly output} 
    20602025 
    2061 %------------------------------------------nam_diatmb----------------------------------------------------- 
    2062  
    2063 \nlst{nam_diatmb} 
    2064 %---------------------------------------------------------------------------------------------------------- 
     2026\begin{listing} 
     2027  \nlst{nam_diatmb} 
     2028  \caption{\forcode{&nam_diatmb}} 
     2029  \label{lst:nam_diatmb} 
     2030\end{listing} 
    20652031 
    20662032A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. 
    … …  
    20702036This diagnostic is actived with the logical $ln\_diatmb$. 
    20712037 
    2072 % ----------------------------------------------------------- 
    2073 %     Courant numbers 
    2074 % ----------------------------------------------------------- 
     2038%% ================================================================================================= 
    20752039\subsection{Courant numbers} 
    20762040 
    … …  
    20802044\[ 
    20812045  C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}} 
    2082   % \label{eq:CFL} 
     2046  % \label{eq:DIA_CFL} 
    20832047\] 
    20842048 
    … …  
    20892053Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step. 
    20902054 
    2091 The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. 
     2055The variables can be activated by setting the \np{nn_diacfl}{nn\_diacfl} namelist parameter to 1 in the \nam{ctl}{ctl} namelist. 
    20922056The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. 
    20932057In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of 
    … …  
    20952059At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along 
    20962060with the coordinates of each. 
    2097 The maximum values from the run are also copied to the ocean.output file.  
    2098  
    2099 % ================================================================ 
    2100  
    2101 \biblio 
    2102  
    2103 \pindex 
     2061The maximum values from the run are also copied to the ocean.output file. 
     2062 
     2063\subinc{\input{../../global/epilogue}} 
    21042064 
    21052065\end{document} 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/subfiles/chap_DIU.tex

    r10442 r11830  
    22 
    33\begin{document} 
    4 % ================================================================ 
    5 % Diurnal SST models (DIU) 
    6 % Edited by James While 
    7 % ================================================================ 
     4 
    85\chapter{Diurnal SST Models (DIU)} 
    96\label{chap:DIU} 
    107 
    11 \minitoc 
     8\thispagestyle{plain} 
    129 
     10\chaptertoc 
    1311 
    14 \newpage 
    15 $\ $\newline % force a new line 
     12\paragraph{Changes record} ~\\ 
     13 
     14{\footnotesize 
     15  \begin{tabularx}{\textwidth}{l||X|X} 
     16    Release & Author(s) & Modifications \\ 
     17    \hline 
     18    {\em   4.0} & {\em ...} & {\em ...} \\ 
     19    {\em   3.6} & {\em ...} & {\em ...} \\ 
     20    {\em   3.4} & {\em ...} & {\em ...} \\ 
     21    {\em <=3.4} & {\em ...} & {\em ...} 
     22  \end{tabularx} 
     23} 
     24 
     25\clearpage 
    1626 
    1727Code to produce an estimate of the diurnal warming and cooling of the sea surface skin 
    18 temperature (skin SST) is found in the DIU directory.   
     28temperature (skin SST) is found in the DIU directory. 
    1929The skin temperature can be split into three parts: 
    2030\begin{itemize} 
    21 \item 
    22   A foundation SST which is free from diurnal warming. 
    23 \item 
    24   A warm layer, typically ~3\,m thick, 
     31\item A foundation SST which is free from diurnal warming. 
     32\item A warm layer, typically ~3\,m thick, 
    2533  where heating from solar radiation can cause a warm stably stratified layer during the daytime 
    26 \item 
    27   A cool skin, a thin layer, approximately ~1\, mm thick, 
     34\item A cool skin, a thin layer, approximately ~1\, mm thick, 
    2835  where long wave cooling is dominant and cools the immediate ocean surface. 
    2936\end{itemize} 
    3037 
    3138Models are provided for both the warm layer, \mdl{diurnal\_bulk}, and the cool skin, \mdl{cool\_skin}. 
    32 Foundation SST is not considered as it can be obtained either from the main NEMO model 
    33 (\ie from the temperature of the top few model levels) or from some other source.   
     39Foundation SST is not considered as it can be obtained either from the main \NEMO\ model 
     40(\ie\ from the temperature of the top few model levels) or from some other source. 
    3441It must be noted that both the cool skin and warm layer models produce estimates of the change in temperature 
    35 ($\Delta T_{\rm{cs}}$ and $\Delta T_{\rm{wl}}$) and 
     42($\Delta T_{\mathrm{cs}}$ and $\Delta T_{\mathrm{wl}}$) and 
    3643both must be added to a foundation SST to obtain the true skin temperature. 
    3744 
    38 Both the cool skin and warm layer models are controlled through the namelist \ngn{namdiu}: 
     45Both the cool skin and warm layer models are controlled through the namelist \nam{diu}{diu}: 
    3946 
    40 \nlst{namdiu} 
     47\begin{listing} 
     48  \nlst{namdiu} 
     49  \caption{\forcode{&namdiu}} 
     50  \label{lst:namdiu} 
     51\end{listing} 
     52 
    4153This namelist contains only two variables: 
    4254\begin{description} 
    43 \item[\np{ln\_diurnal}] 
    44   A logical switch for turning on/off both the cool skin and warm layer. 
    45 \item[\np{ln\_diurnal\_only}] 
    46   A logical switch which if \forcode{.true.} will run the diurnal model without the other dynamical parts of NEMO. 
    47   \np{ln\_diurnal\_only} must be \forcode{.false.} if \np{ln\_diurnal} is \forcode{.false.}. 
     55\item [{\np{ln_diurnal}{ln\_diurnal}}] A logical switch for turning on/off both the cool skin and warm layer. 
     56\item [{\np{ln_diurnal_only}{ln\_diurnal\_only}}] A logical switch which if \forcode{.true.} will run the diurnal model without the other dynamical parts of \NEMO. 
     57  \np{ln_diurnal_only}{ln\_diurnal\_only} must be \forcode{.false.} if \np{ln_diurnal}{ln\_diurnal} is \forcode{.false.}. 
    4858\end{description} 
    4959 
    … …  
    5363Initialisation is through the restart file. 
    5464Specifically the code will expect the presence of the 2-D variable ``Dsst'' to initialise the warm layer. 
    55 The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable.   
     65The cool skin model, which is determined purely by the instantaneous fluxes, has no initialisation variable. 
    5666 
    57 %=============================================================== 
     67%% ================================================================================================= 
    5868\section{Warm layer model} 
    59 \label{sec:warm_layer_sec} 
    60 %=============================================================== 
     69\label{sec:DIU_warm_layer_sec} 
    6170 
    62 The warm layer is calculated using the model of \citet{Takaya_al_JGR10} (TAKAYA10 model hereafter). 
     71The warm layer is calculated using the model of \citet{takaya.bidlot.ea_JGR10} (TAKAYA10 model hereafter). 
    6372This is a simple flux based model that is defined by the equations 
    6473\begin{align} 
    65 \frac{\partial{\Delta T_{\rm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 
     74\frac{\partial{\Delta T_{\mathrm{wl}}}}{\partial{t}}&=&\frac{Q(\nu+1)}{D_T\rho_w c_p 
    6675\nu}-\frac{(\nu+1)ku^*_{w}f(L_a)\Delta T}{D_T\Phi\!\left(\frac{D_T}{L}\right)} \mbox{,} 
    67 \label{eq:ecmwf1} \\ 
    68 L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:ecmwf2} 
     76\label{eq:DIU_ecmwf1} \\ 
     77L&=&\frac{\rho_w c_p u^{*^3}_{w}}{\kappa g \alpha_w Q }\mbox{,}\label{eq:DIU_ecmwf2} 
    6978\end{align} 
    70 where $\Delta T_{\rm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. 
    71 In equation (\autoref{eq:ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, 
     79where $\Delta T_{\mathrm{wl}}$ is the temperature difference between the top of the warm layer and the depth $D_T=3$\,m at which there is assumed to be no diurnal signal. 
     80In equation (\autoref{eq:DIU_ecmwf1}) $\alpha_w=2\times10^{-4}$ is the thermal expansion coefficient of water, 
    7281$\kappa=0.4$ is von K\'{a}rm\'{a}n's constant, $c_p$ is the heat capacity at constant pressure of sea water, 
    7382$\rho_w$ is the water density, and $L$ is the Monin-Obukhov length. 
    7483The tunable variable $\nu$ is a shape parameter that defines the expected subskin temperature profile via 
    75 $T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\rm{wl}}$, 
     84$T(z) = T(0) - \left( \frac{z}{D_T} \right)^\nu \Delta T_{\mathrm{wl}}$, 
    7685where $T$ is the absolute temperature and $z\le D_T$ is the depth below the top of the warm layer. 
    7786The influence of wind on TAKAYA10 comes through the magnitude of the friction velocity of the water $u^*_{w}$, 
    … …  
    7988the relationship $u^*_{w} = u_{10}\sqrt{\frac{C_d\rho_a}{\rho_w}}$, where $C_d$ is the drag coefficient, 
    8089and $\rho_a$ is the density of air. 
    81 The symbol $Q$ in equation (\autoref{eq:ecmwf1}) is the instantaneous total thermal energy flux into 
     90The symbol $Q$ in equation (\autoref{eq:DIU_ecmwf1}) is the instantaneous total thermal energy flux into 
    8291the diurnal layer, \ie 
    8392\[ 
    84   Q = Q_{\rm{sol}} + Q_{\rm{lw}} + Q_{\rm{h}}\mbox{,} 
    85   % \label{eq:e_flux_eqn} 
     93  Q = Q_{\mathrm{sol}} + Q_{\mathrm{lw}} + Q_{\mathrm{h}}\mbox{,} 
     94  % \label{eq:DIU_e_flux_eqn} 
    8695\] 
    87 where $Q_{\rm{h}}$ is the sensible and latent heat flux, $Q_{\rm{lw}}$ is the long wave flux, 
    88 and $Q_{\rm{sol}}$ is the solar flux absorbed within the diurnal warm layer. 
    89 For $Q_{\rm{sol}}$ the 9 term representation of \citet{Gentemann_al_JGR09} is used. 
    90 In equation \autoref{eq:ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, 
     96where $Q_{\mathrm{h}}$ is the sensible and latent heat flux, $Q_{\mathrm{lw}}$ is the long wave flux, 
     97and $Q_{\mathrm{sol}}$ is the solar flux absorbed within the diurnal warm layer. 
     98For $Q_{\mathrm{sol}}$ the 9 term representation of \citet{gentemann.minnett.ea_JGR09} is used. 
     99In equation \autoref{eq:DIU_ecmwf1} the function $f(L_a)=\max(1,L_a^{\frac{2}{3}})$, 
    91100where $L_a=0.3$\footnote{ 
    92101  This is a global average value, more accurately $L_a$ could be computed as $L_a=(u^*_{w}/u_s)^{\frac{1}{2}}$, 
    … …  
    991084\zeta^2}{1+3\zeta+0.25\zeta^2} &(\zeta \ge 0) \\ 
    100109                                    (1 - 16\zeta)^{-\frac{1}{2}} & (\zeta < 0) \mbox{,} 
    101                                     \end{array} \right. \label{eq:stab_func_eqn} 
     110                                    \end{array} \right. \label{eq:DIU_stab_func_eqn} 
    102111\end{equation} 
    103 where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:stab_func_eqn}), 
    104 and thus of (\autoref{eq:ecmwf1}), is discontinuous at $\zeta=0$ (\ie $Q\rightarrow0$ in 
    105 equation (\autoref{eq:ecmwf2})). 
     112where $\zeta=\frac{D_T}{L}$.  It is clear that the first derivative of (\autoref{eq:DIU_stab_func_eqn}), 
     113and thus of (\autoref{eq:DIU_ecmwf1}), is discontinuous at $\zeta=0$ (\ie\ $Q\rightarrow0$ in 
     114equation (\autoref{eq:DIU_ecmwf2})). 
    106115 
    107 The two terms on the right hand side of (\autoref{eq:ecmwf1}) represent different processes. 
     116The two terms on the right hand side of (\autoref{eq:DIU_ecmwf1}) represent different processes. 
    108117The first term is simply the diabatic heating or cooling of the diurnal warm layer due to 
    109118thermal energy fluxes into and out of the layer. 
    … …  
    111120In practice the second term acts as a relaxation on the temperature. 
    112121 
    113 %=============================================================== 
     122%% ================================================================================================= 
     123\section{Cool skin model} 
     124\label{sec:DIU_cool_skin_sec} 
    114125 
    115 \section{Cool skin model} 
    116 \label{sec:cool_skin_sec} 
    117  
    118 %=============================================================== 
    119  
    120 The cool skin is modelled using the framework of \citet{Saunders_JAS82} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$. 
    121 As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\rm{cs}}$ becomes 
     126The cool skin is modelled using the framework of \citet{saunders_JAS67} who used a formulation of the near surface temperature difference based upon the heat flux and the friction velocity $u^*_{w}$. 
     127As the cool skin is so thin (~1\,mm) we ignore the solar flux component to the heat flux and the Saunders equation for the cool skin temperature difference $\Delta T_{\mathrm{cs}}$ becomes 
    122128\[ 
    123   % \label{eq:sunders_eqn} 
    124   \Delta T_{\rm{cs}}=\frac{Q_{\rm{ns}}\delta}{k_t} \mbox{,} 
     129  % \label{eq:DIU_sunders_eqn} 
     130  \Delta T_{\mathrm{cs}}=\frac{Q_{\mathrm{ns}}\delta}{k_t} \mbox{,} 
    125131\] 
    126 where $Q_{\rm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and 
     132where $Q_{\mathrm{ns}}$ is the, usually negative, non-solar heat flux into the ocean and 
    127133$k_t$ is the thermal conductivity of sea water. 
    128134$\delta$ is the thickness of the skin layer and is given by 
    129135\begin{equation} 
    130 \label{eq:sunders_thick_eqn} 
     136\label{eq:DIU_sunders_thick_eqn} 
    131137\delta=\frac{\lambda \mu}{u^*_{w}} \mbox{,} 
    132138\end{equation} 
    133139where $\mu$ is the kinematic viscosity of sea water and $\lambda$ is a constant of proportionality which 
    134 \citet{Saunders_JAS82} suggested varied between 5 and 10. 
     140\citet{saunders_JAS67} suggested varied between 5 and 10. 
    135141 
    136 The value of $\lambda$ used in equation (\autoref{eq:sunders_thick_eqn}) is that of \citet{Artale_al_JGR02}, 
    137 which is shown in \citet{Tu_Tsuang_GRL05} to outperform a number of other parametrisations at 
     142The value of $\lambda$ used in equation (\autoref{eq:DIU_sunders_thick_eqn}) is that of \citet{artale.iudicone.ea_JGR02}, 
     143which is shown in \citet{tu.tsuang_GRL05} to outperform a number of other parametrisations at 
    138144both low and high wind speeds. 
    139145Specifically, 
    140146\[ 
    141   % \label{eq:artale_lambda_eqn} 
     147  % \label{eq:DIU_artale_lambda_eqn} 
    142148  \lambda = \frac{ 8.64\times10^4 u^*_{w} k_t }{ \rho c_p h \mu \gamma }\mbox{,} 
    143149\] 
    … …  
    145151$\gamma$ is a dimensionless function of wind speed $u$: 
    146152\[ 
    147   % \label{eq:artale_gamma_eqn} 
     153  % \label{eq:DIU_artale_gamma_eqn} 
    148154  \gamma = 
    149155  \begin{cases} 
    … …  
    154160\] 
    155161 
    156 \biblio 
    157  
    158 \pindex 
     162\subinc{\input{../../global/epilogue}} 
    159163 
    160164\end{document} 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/subfiles/chap_DOM.tex

    r10502 r11830  
    22 
    33\begin{document} 
    4 % ================================================================ 
    5 % Chapter 2 ——— Space and Time Domain (DOM) 
    6 % ================================================================ 
     4 
    75\chapter{Space Domain (DOM)} 
    86\label{chap:DOM} 
    97 
    10 \minitoc 
    11  
    12 % Missing things: 
    13 %  - istate: description of the initial state   ==> this has to be put elsewhere.. 
    14 %                  perhaps in MISC ?  By the way the initialisation of T S and dynamics  
    15 %                  should be put outside of DOM routine (better with TRC staff and off-line 
    16 %                  tracers) 
    17 %  -geo2ocean:  how to switch from geographic to mesh coordinate 
    18 %     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 
    19  
    20 \newpage 
    21  
    22 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, 
    23 we need to choose a discretization on a grid, and numerical algorithms. 
     8% Missing things 
     9% -    istate: description of the initial state   ==> this has to be put elsewhere.. 
     10%              perhaps in MISC ?  By the way the initialisation of T S and dynamics 
     11%              should be put outside of DOM routine (better with TRC staff and off-line 
     12%              tracers) 
     13% - geo2ocean: how to switch from geographic to mesh coordinate 
     14% -    domclo: closed sea and lakes.... 
     15%              management of closea sea area: specific to global cfg, both forced and coupled 
     16 
     17\thispagestyle{plain} 
     18 
     19\chaptertoc 
     20 
     21\paragraph{Changes record} ~\\ 
     22 
     23{\footnotesize 
     24  \begin{tabularx}{0.8\textwidth}{l||X|X} 
     25    Release                                                                                 & 
     26    Author(s)                                                                               & 
     27    Modifications                                                                           \\ 
     28    \hline 
     29    {\em 4.0                                                                              } & 
     30    {\em Simon M\"{u}ller \& Andrew Coward \newline \newline 
     31      Simona Flavoni and Tim Graham                                                       } & 
     32    {\em Compatibility changes: many options moved to external domain configuration tools 
     33      (see \autoref{apdx:DOMCFG}). \newline 
     34      Updates                                                                             } \\ 
     35    {\em 3.6                                                                              } & 
     36    {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec         } & 
     37    {\em Updates                                                                          } \\ 
     38    {\em $\leq$ 3.4                                                                       } & 
     39    {\em Gurvan Madec and S\'{e}bastien Masson                                            } & 
     40    {\em First version                                                                    } 
     41  \end{tabularx} 
     42} 
     43 
     44\clearpage 
     45 
     46Having defined the continuous equations in \autoref{chap:MB} and 
     47chosen a time discretisation \autoref{chap:TD}, 
     48we need to choose a grid for spatial discretisation and related numerical algorithms. 
    2449In the present chapter, we provide a general description of the staggered grid used in \NEMO, 
    25 and other information relevant to the main directory routines as well as the DOM (DOMain) directory. 
    26  
    27 % ================================================================ 
    28 % Fundamentals of the Discretisation 
    29 % ================================================================ 
     50and other relevant information about the DOM (DOMain) source code modules. 
     51 
     52%% ================================================================================================= 
    3053\section{Fundamentals of the discretisation} 
    3154\label{sec:DOM_basics} 
    3255 
    33 % ------------------------------------------------------------------------------------------------------------- 
    34 %        Arrangement of Variables  
    35 % ------------------------------------------------------------------------------------------------------------- 
     56%% ================================================================================================= 
    3657\subsection{Arrangement of variables} 
    3758\label{subsec:DOM_cell} 
    3859 
    39 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    40 \begin{figure}[!tb] 
    41   \begin{center} 
    42     \includegraphics[]{Fig_cell} 
    43     \caption{ 
    44       \protect\label{fig:cell} 
    45       Arrangement of variables. 
    46       $t$ indicates scalar points where temperature, salinity, density, pressure and 
    47       horizontal divergence are defined. 
    48       $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and 
    49       planetary vorticities are defined. 
    50     } 
    51   \end{center} 
     60\begin{figure} 
     61  \centering 
     62  \includegraphics[width=0.33\textwidth]{DOM_cell} 
     63  \caption[Arrangement of variables in the unit cell of space domain]{ 
     64    Arrangement of variables in the unit cell of space domain. 
     65    $t$ indicates scalar points where 
     66    temperature, salinity, density, pressure and horizontal divergence are defined. 
     67    $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where 
     68    both relative and planetary vorticities are defined.} 
     69  \label{fig:DOM_cell} 
    5270\end{figure} 
    53 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    54  
    55 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 
    56 centred second-order finite difference approximation. 
    57 Special attention has been given to the homogeneity of the solution in the three space directions. 
     71 
     72The numerical techniques used to solve the Primitive Equations in this model are based on 
     73the traditional, centred second-order finite difference approximation. 
     74Special attention has been given to the homogeneity of the solution in the three spatial directions. 
    5875The arrangement of variables is the same in all directions. 
    59 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 
    60 the centre of each face of the cells (\autoref{fig:cell}). 
    61 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 
    62 \citep{Mesinger_Arakawa_Bk76}. 
    63 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 
    64 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 
    65  
    66 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 
    67 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    68 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. 
    69 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 
    70 the grid-point where the scale factors are defined. 
    71 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 
     76It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with 
     77vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). 
     78This is the generalisation to three dimensions of the well-known ``C'' grid in 
     79Arakawa's classification \citep{mesinger.arakawa_bk76}. 
     80The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each 
     81vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying 
     82the $\zeta$ and $f$-points. 
     83 
     84The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by 
     85the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
     86The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on 
     87\autoref{tab:DOM_cell}. 
     88In all the following, 
     89subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where 
     90the scale factors are defined. 
     91Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 
    7292As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
    73 $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. 
    74 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 
    75 while the scale factors are chosen equal to their local analytical value. 
     93$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 
     94Discrete partial derivatives are formulated by 
     95the traditional, centred second order finite difference approximation while 
     96the scale factors are chosen equal to their local analytical value. 
    7697An important point here is that the partial derivative of the scale factors must be evaluated by 
    7798centred finite difference approximation, not from their analytical expression. 
    78 This preserves the symmetry of the discrete set of equations and therefore satisfies many of 
    79 the continuous properties (see \autoref{apdx:C}). 
     99This preserves the symmetry of the discrete set of equations and 
     100therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). 
    80101A similar, related remark can be made about the domain size: 
    81 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors 
    82 (see \autoref{eq:DOM_bar} in the next section). 
    83  
    84 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    85 \begin{table}[!tb] 
    86   \begin{center} 
    87     \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 
    88       \hline 
    89       T  & $i      $ & $j      $ & $k      $ \\ 
    90       \hline 
    91       u  & $i + 1/2$ & $j      $ & $k      $ \\ 
    92       \hline 
    93       v  & $i      $ & $j + 1/2$ & $k      $ \\ 
    94       \hline 
    95       w  & $i      $ & $j      $ & $k + 1/2$ \\ 
    96       \hline 
    97       f  & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
    98       \hline 
    99       uw & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
    100       \hline 
    101       vw & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
    102       \hline 
    103       fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
    104       \hline 
    105     \end{tabular} 
    106     \caption{ 
    107       \protect\label{tab:cell} 
    108       Location of grid-points as a function of integer or integer and a half value of the column, line or level. 
    109       This indexing is only used for the writing of the semi -discrete equation. 
    110       In the code, the indexing uses integer values only and has a reverse direction in the vertical 
    111       (see \autoref{subsec:DOM_Num_Index}) 
    112     } 
    113   \end{center} 
     102when needed, an area, volume, or the total ocean depth must be evaluated as 
     103the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). 
     104 
     105\begin{table} 
     106  \centering 
     107  \begin{tabular}{|l|l|l|l|} 
     108    \hline 
     109    t   & $i      $ & $j      $ & $k      $ \\ 
     110    \hline 
     111    u   & $i + 1/2$ & $j      $ & $k      $ \\ 
     112    \hline 
     113    v   & $i      $ & $j + 1/2$ & $k      $ \\ 
     114    \hline 
     115    w   & $i      $ & $j      $ & $k + 1/2$ \\ 
     116    \hline 
     117    f   & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
     118    \hline 
     119    uw  & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
     120    \hline 
     121    vw  & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
     122    \hline 
     123    fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
     124    \hline 
     125  \end{tabular} 
     126  \caption[Location of grid-points]{ 
     127    Location of grid-points as a function of integer or 
     128    integer and a half value of the column, line or level. 
     129    This indexing is only used for the writing of the semi-discrete equations. 
     130    In the code, the indexing uses integer values only and 
     131    is positive downwards in the vertical with $k=1$ at the surface. 
     132    (see \autoref{subsec:DOM_Num_Index})} 
     133  \label{tab:DOM_cell} 
    114134\end{table} 
    115 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    116  
    117 % ------------------------------------------------------------------------------------------------------------- 
    118 %        Vector Invariant Formulation  
    119 % ------------------------------------------------------------------------------------------------------------- 
     135 
     136Note that the definition of the scale factors 
     137(\ie\ as the analytical first derivative of the transformation that 
     138results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 
     139is specific to the \NEMO\ model \citep{marti.madec.ea_JGR92}. 
     140As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, 
     141whereas many other models on a C grid choose to define such a scale factor as 
     142the distance between the $u$-points on each side of the $t$-point. 
     143Relying on an analytical transformation has two advantages: 
     144firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
     145since they are first introduced in the continuous equations; 
     146secondly, analytical transformations encourage good practice by 
     147the definition of smoothly varying grids 
     148(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) 
     149\citep{treguier.dukowicz.ea_JGR96}. 
     150An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 
     151\begin{figure} 
     152  \centering 
     153  \includegraphics[width=0.5\textwidth]{DOM_zgr_e3} 
     154  \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ 
     155    Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
     156    and (b) analytically derived grid-point position and scale factors. 
     157    For both grids here, the same $w$-point depth has been chosen but 
     158    in (a) the $t$-points are set half way between $w$-points while 
     159    in (b) they are defined from an analytical function: 
     160    $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 
     161    Note the resulting difference between the value of the grid-size $\Delta_k$ and 
     162    those of the scale factor $e_k$.} 
     163  \label{fig:DOM_zgr_e3} 
     164\end{figure} 
     165 
     166%% ================================================================================================= 
    120167\subsection{Discrete operators} 
    121168\label{subsec:DOM_operators} 
    122169 
    123 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at 
    124 the midpoint between them are: 
     170Given the values of a variable $q$ at adjacent points, 
     171the differencing and averaging operators at the midpoint between them are: 
    125172\begin{alignat*}{2} 
    126   % \label{eq:di_mi} 
     173  % \label{eq:DOM_di_mi} 
    127174  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\ 
    128175  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 
    … …  
    130177 
    131178Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 
    132 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 
    133 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 
    134 its Laplacian is defined at $t$-point. 
     179Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, 
     180the gradient of a variable $q$ defined at a $t$-point has 
     181its three components defined at $u$-, $v$- and $w$-points while 
     182its Laplacian is defined at the $t$-point. 
    135183These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
    136 \[ 
     184\begin{gather*} 
    137185  % \label{eq:DOM_grad} 
    138186  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i 
    139187                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j 
    140                   + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 
    141 \] 
    142 \begin{multline*} 
     188                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\ 
    143189  % \label{eq:DOM_lap} 
    144190  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    145191                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) 
    146                             + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ 
     192                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] 
    147193                  + \frac{1}{e_{3t}} 
    148194                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 
    149 \end{multline*} 
    150  
    151 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
    152 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 
     195\end{gather*} 
     196 
     197Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, 
     198a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has 
     199its three curl components defined at $vw$-, $uw$, and $f$-points, and 
    153200its divergence defined at $t$-points: 
    154 \begin{multline} 
     201\begin{multline*} 
    155202% \label{eq:DOM_curl} 
    156203  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} 
    … …  
    163210                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) 
    164211                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 
    165 \end{multline} 
    166 \begin{equation} 
     212\end{multline*} 
     213\[ 
    167214% \label{eq:DOM_div} 
    168215  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    169216                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 
    170217                              + \frac{1}{e_{3t}} \delta_k (a_3) 
    171 \end{equation} 
    172  
    173 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which 
    174 is a masked field (i.e. equal to zero inside solid area): 
     218\] 
     219 
     220The vertical average over the whole water column is denoted by an overbar and 
     221is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 
    175222\begin{equation} 
    176223  \label{eq:DOM_bar} 
    … …  
    178225\end{equation} 
    179226where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 
    180 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over 
    181 all grid points of the same type in the direction indicated by the subscript (here $k$). 
     227$k^b$ and $k^o$ are the bottom and surface $k$-indices, 
     228and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in 
     229the direction indicated by the subscript (here $k$). 
    182230 
    183231In continuous form, the following properties are satisfied: 
    … …  
    189237\end{gather} 
    190238 
    191 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as 
    192 the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 
     239It is straightforward to demonstrate that these properties are verified locally in discrete form as 
     240soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 
    193241vector points $(u,v,w)$. 
    194242 
    195 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. 
    196 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
    197 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 
    198 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
    199 \begin{alignat}{4} 
     243Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 
     244It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and 
     245$\delta_k$) are skew-symmetric linear operators, 
     246and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and 
     247$\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
     248\begin{alignat}{5} 
    200249  \label{eq:DOM_di_adj} 
    201250  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\ 
    … …  
    204253\end{alignat} 
    205254 
    206 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and  
     255In other words, 
     256the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 
    207257$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 
    208 These two properties will be used extensively in the \autoref{apdx:C} to 
     258These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to 
    209259demonstrate integral conservative properties of the discrete formulation chosen. 
    210260 
    211 % ------------------------------------------------------------------------------------------------------------- 
    212 %        Numerical Indexing  
    213 % ------------------------------------------------------------------------------------------------------------- 
     261%% ================================================================================================= 
    214262\subsection{Numerical indexing} 
    215263\label{subsec:DOM_Num_Index} 
    216264 
    217 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    218 \begin{figure}[!tb] 
    219   \begin{center} 
    220     \includegraphics[]{Fig_index_hor} 
    221     \caption{ 
    222       \protect\label{fig:index_hor} 
    223       Horizontal integer indexing used in the \fortran code. 
    224       The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices 
    225     } 
    226   \end{center} 
     265\begin{figure} 
     266  \centering 
     267  \includegraphics[width=0.33\textwidth]{DOM_index_hor} 
     268  \caption[Horizontal integer indexing]{ 
     269    Horizontal integer indexing used in the \fortran\ code. 
     270    The dashed area indicates the cell in which 
     271    variables contained in arrays have the same $i$- and $j$-indices} 
     272  \label{fig:DOM_index_hor} 
    227273\end{figure} 
    228 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    229  
    230 The array representation used in the \fortran code requires an integer indexing while 
    231 the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 
    232 integer values for $t$-points and both integer and integer and a half values for all the other points. 
    233 Therefore a specific integer indexing must be defined for points other than $t$-points 
    234 (\ie velocity and vorticity grid-points). 
    235 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. 
    236  
    237 % ----------------------------------- 
    238 %        Horizontal Indexing  
    239 % ----------------------------------- 
     274 
     275The array representation used in the \fortran\ code requires an integer indexing. 
     276However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with 
     277the use of integer values for $t$-points only while 
     278all the other points involve integer and a half values. 
     279Therefore, a specific integer indexing has been defined for points other than $t$-points 
     280(\ie\ velocity and vorticity grid-points). 
     281Furthermore, the direction of the vertical indexing has been reversed and 
     282the surface level set at $k = 1$. 
     283 
     284%% ================================================================================================= 
    240285\subsubsection{Horizontal indexing} 
    241286\label{subsec:DOM_Num_Index_hor} 
    242287 
    243 The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}. 
     288The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}. 
    244289For an increasing $i$ index ($j$ index), 
    245290the $t$-point and the eastward $u$-point (northward $v$-point) have the same index 
    246 (see the dashed area in \autoref{fig:index_hor}). 
    247 A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. 
    248  
    249 % ----------------------------------- 
    250 %        Vertical indexing  
    251 % ----------------------------------- 
     291(see the dashed area in \autoref{fig:DOM_index_hor}). 
     292A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices. 
     293 
     294%% ================================================================================================= 
    252295\subsubsection{Vertical indexing} 
    253296\label{subsec:DOM_Num_Index_vertical} 
    254297 
    255 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in 
    256 the \fortran code compared to the indexing used in the semi -discrete equations and 
    257 given in \autoref{subsec:DOM_cell}. 
    258 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below 
    259 (\autoref{fig:index_vert}). 
    260 The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 
    261 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 
    262 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, 
    263 in opposition to what is done in the horizontal plane where 
    264 it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that 
    265 have the same $i$ or $j$ index 
    266 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 
     298In the vertical, the chosen indexing requires special attention since 
     299the direction of the $k$-axis in the \fortran\ code is the reverse of 
     300that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 
     301The sea surface corresponds to the $w$-level $k = 1$, 
     302which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). 
     303The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 
     304the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 
     305Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 
     306(\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 
     307in contrast to the indexing on the horizontal plane where 
     308the $t$-point has the same index as the nearest velocity points in 
     309the positive direction of the respective horizontal axis index 
     310(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 
    267311Since the scale factors are chosen to be strictly positive, 
    268 a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of 
    269 the discrete equations given in this documentation. 
    270  
    271 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    272 \begin{figure}[!pt] 
    273   \begin{center} 
    274     \includegraphics[]{Fig_index_vert} 
    275     \caption{ 
    276       \protect\label{fig:index_vert} 
    277       Vertical integer indexing used in the \fortran code. 
    278       Note that the $k$-axis is orientated downward. 
    279       The dashed area indicates the cell in which variables contained in arrays have the same $k$-index. 
    280     } 
    281   \end{center} 
     312a \textit{minus sign} is included in the \fortran\ implementations of 
     313\textit{all the vertical derivatives} of the discrete equations given in this manual in order to 
     314accommodate the opposing vertical index directions in implementation and documentation. 
     315 
     316\begin{figure} 
     317  \centering 
     318  \includegraphics[width=0.33\textwidth]{DOM_index_vert} 
     319  \caption[Vertical integer indexing]{ 
     320    Vertical integer indexing used in the \fortran\ code. 
     321    Note that the $k$-axis is oriented downward. 
     322    The dashed area indicates the cell in which 
     323    variables contained in arrays have a common $k$-index.} 
     324  \label{fig:DOM_index_vert} 
    282325\end{figure} 
    283 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    284  
    285 % ----------------------------------- 
    286 %        Domain Size 
    287 % ----------------------------------- 
    288 \subsubsection{Domain size} 
     326 
     327%% ================================================================================================= 
     328\section{Spatial domain configuration} 
     329\label{subsec:DOM_config} 
     330 
     331Two typical methods are available to specify the spatial domain configuration; 
     332they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in 
     333namelist \nam{cfg}{cfg}. 
     334 
     335If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, 
     336the domain-specific parameters and fields are read from a NetCDF input file, 
     337whose name (without its .nc suffix) can be specified as 
     338the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 
     339 
     340If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, 
     341the domain-specific parameters and fields can be provided (\eg\ analytically computed) by 
     342subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
     343These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, 
     344and default versions that configure the spatial domain for the GYRE reference configuration are 
     345present in the \path{./src/OCE/USR} directory. 
     346 
     347In version 4.0 there are no longer any options for reading complex bathymetries and 
     348performing a vertical discretisation at run-time. 
     349Whilst it is occasionally convenient to have a common bathymetry file and, for example, 
     350to run similar models with and without partial bottom boxes and/or sigma-coordinates, 
     351supporting such choices leads to overly complex code. 
     352Worse still is the difficulty of ensuring the model configurations intended to be identical are 
     353indeed so when the model domain itself can be altered by runtime selections. 
     354The code previously used to perform vertical discretisation has been incorporated into 
     355an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 
     356 
     357The next subsections summarise the parameter and fields related to 
     358the configuration of the whole model domain. 
     359These represent the minimum information that must be provided either via 
     360the \np{cn_domcfg}{cn\_domcfg} file or 
     361set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 
     362The requirements are presented in three sections: 
     363the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), 
     364and the vertical grid (\autoref{subsec:DOM_zgr}). 
     365 
     366%% ================================================================================================= 
     367\subsection{Domain size} 
    289368\label{subsec:DOM_size} 
    290369 
    291 The total size of the computational domain is set by the parameters \np{jpiglo}, 
    292 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 
    293 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 
    294 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 
    295 see \autoref{sec:LBC_mpp}). 
    296  
    297 % ================================================================ 
    298 % Domain: List of fields needed 
    299 % ================================================================ 
    300 \section{Needed fields} 
    301 \label{sec:DOM_fields} 
    302 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 
    303 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    304 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 
    305 The associated scale factors are defined using the analytical first derivative of the transformation 
    306 \autoref{eq:scale_factors}. 
    307 Necessary fields for configuration definition are: 
     370The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and 
     371\jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. 
     372Note, that the variables \texttt{jpi} and \texttt{jpj} refer to 
     373the size of each processor subdomain when the code is run in parallel using domain decomposition 
     374(\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). 
     375 
     376The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, 
     377and the nominal resolution through parameter \np{nn_cfg}{nn\_cfg} 
     378(unless in the input file both of variables \texttt{ORCA} and \texttt{ORCA\_index} are present, 
     379in which case \np{cn_cfg}{cn\_cfg} and \np{nn_cfg}{nn\_cfg} are set from these values accordingly). 
     380 
     381The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. 
     382See \autoref{sec:LBC_jperio} for details on the available options and 
     383the corresponding values for \jp{jperio}. 
     384 
     385%% ================================================================================================= 
     386\subsection[Horizontal grid mesh (\textit{domhgr.F90}]{Horizontal grid mesh (\protect\mdl{domhgr})} 
     387\label{subsec:DOM_hgr} 
     388 
     389%% ================================================================================================= 
     390\subsubsection{Required fields} 
     391\label{sec:DOM_hgr_fields} 
     392 
     393The explicit specification of a range of mesh-related fields are required for 
     394the definition of a configuration. 
     395These include: 
     396 
     397\begin{clines} 
     398int    jpiglo, jpjglo, jpkglo     /* global domain sizes                                    */ 
     399int    jperio                     /* lateral global domain b.c.                             */ 
     400double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ 
     401double gphit, gphiu, gphiv, gphif /* geographic latitude                                    */ 
     402double e1t, e1u, e1v, e1f         /* horizontal scale factors                               */ 
     403double e2t, e2u, e2v, e2f         /* horizontal scale factors                               */ 
     404\end{clines} 
     405 
     406The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to 
     407the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, 
     408evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position. 
     409The calculation of the values of the horizontal scale factor arrays in general additionally involves 
     410partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, 
     411evaluated for the same arguments as $\lambda$ and $\varphi$. 
     412 
     413%% ================================================================================================= 
     414\subsubsection{Optional fields} 
     415 
     416\begin{clines} 
     417                        /* Optional:                                                 */ 
     418int    ORCA, ORCA_index /* configuration name, configuration resolution              */ 
     419double e1e2u, e1e2v     /* U and V surfaces (if grid size reduction in some straits) */ 
     420double ff_f, ff_t       /* Coriolis parameter (if not on the sphere)                 */ 
     421\end{clines} 
     422 
     423\NEMO\ can support the local reduction of key strait widths by 
     424altering individual values of e2u or e1v at the appropriate locations. 
     425This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points 
     426(see \autoref{sec:MISC_strait} for illustrated examples). 
     427The key is to reduce the faces of $T$-cell 
     428(\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 
     429not the volume of the cells. 
     430Doing otherwise can lead to numerical instability issues. 
     431In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but 
     432in cases where a gridsize reduction is required, 
     433the unaltered surface areas at $u$ and $v$ grid points 
     434(\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. 
     435If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and 
     436the internal computation is suppressed. 
     437Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should 
     438set the surface-area computation flag: 
     439\texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 
     440 
     441\smallskip 
     442Similar logic applies to the other optional fields: 
     443\texttt{ff\_f} and \texttt{ff\_t} which can be used to 
     444provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. 
     445If present these fields will be read and used and 
     446the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 
     447Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should 
     448set the Coriolis computation flag: 
     449\texttt{iff} to a non-zero value to suppress their re-computation. 
     450 
     451Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to 
     452those of $t$ points, thus no specific arrays are defined at $w$ points. 
     453 
     454%% ================================================================================================= 
     455\subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} 
     456\label{subsec:DOM_zgr} 
     457 
     458\begin{listing} 
     459  \nlst{namdom} 
     460  \caption{\forcode{&namdom}} 
     461  \label{lst:namdom} 
     462\end{listing} 
     463 
     464In the vertical, the model mesh is determined by four things: 
     465\begin{enumerate} 
     466\item the bathymetry given in meters; 
     467\item the number of levels of the model (\jp{jpk}); 
     468\item the analytical transformation $z(i,j,k)$ and the vertical scale factors 
     469  (derivatives of the transformation); and 
     470\item the masking system, 
     471  \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. 
     472\end{enumerate} 
     473 
     474\begin{figure} 
     475  \centering 
     476  \includegraphics[width=0.5\textwidth]{DOM_z_zps_s_sps} 
     477  \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ 
     478    The ocean bottom as seen by the model: 
     479    \begin{enumerate*}[label=(\textit{\alph*})] 
     480    \item $z$-coordinate with full step, 
     481    \item $z$-coordinate with partial step, 
     482    \item $s$-coordinate: terrain following representation, 
     483    \item hybrid $s-z$ coordinate, 
     484    \item hybrid $s-z$ coordinate with partial step, and 
     485    \item same as (e) but in the non-linear free surface 
     486      (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 
     487  \end{enumerate*} 
     488  Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 
     489  \label{fig:DOM_z_zps_s_sps} 
     490\end{figure} 
     491 
     492The choice of a vertical coordinate is made when setting up the configuration; 
     493it is not intended to be an option which can be changed in the middle of an experiment. 
     494The one exception to this statement being the choice of linear or non-linear free surface. 
     495In v4.0 the linear free surface option is implemented as 
     496a special case of the non-linear free surface. 
     497This is computationally wasteful since it uses the structures for time-varying 3D metrics 
     498for fields that (in the linear free surface case) are fixed. 
     499However, the linear free-surface is rarely used and 
     500implementing it this way means a single configuration file can support both options. 
     501 
     502By default a non-linear free surface is used 
     503(\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 
     504the coordinate follow the time-variation of the free surface so that 
     505the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 
     506When a linear free surface is assumed 
     507(\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 
     508the vertical coordinates are fixed in time, but 
     509the seawater can move up and down across the $z_0$ surface 
     510(in other words, the top of the ocean in not a rigid lid). 
     511 
     512Note that settings: 
     513\np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} 
     514mentioned in the following sections appear to be namelist options but 
     515they are no longer truly namelist options for \NEMO. 
     516Their value is written to and read from the domain configuration file and 
     517they should be treated as fixed parameters for a particular configuration. 
     518They are namelist options for the \texttt{DOMAINcfg} tool that can be used to 
     519build the configuration file and serve both to provide a record of the choices made whilst 
     520building the configuration and to trigger appropriate code blocks within \NEMO. 
     521These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. 
     522 
     523\medskip 
     524The decision on these choices must be made when the \np{cn_domcfg}{cn\_domcfg} file is constructed. 
     525Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c): 
    308526 
    309527\begin{itemize} 
    310 \item 
    311   Geographic position: 
    312   longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 
    313   latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 
    314   (all respectively at T, U, V and F point) 
    315 \item 
    316   Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 
    317   (at T and F point) 
    318 \item 
    319   Scale factors: 
    320   \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 
    321   \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 
    322   \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 
    323   \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits),  
    324   \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 
     528\item $z$-coordinate with full step bathymetry (\np[=.true.]{ln_zco}{ln\_zco}), 
     529\item $z$-coordinate with partial step ($zps$) bathymetry (\np[=.true.]{ln_zps}{ln\_zps}), 
     530\item Generalized, $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}). 
    325531\end{itemize} 
    326   
    327 These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 
    328 \ngn{namcfg}. 
    329  
    330 \nlst{namcfg} 
    331  
    332 Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 
    333 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 
    334 For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
    335 They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 
    336 they provide the horizontal and vertical mesh. 
    337 % ------------------------------------------------------------------------------------------------------------- 
    338 %        Needed fields  
    339 % ------------------------------------------------------------------------------------------------------------- 
    340 %\subsection{List of needed fields to build DOMAIN} 
    341 %\label{subsec:DOM_fields_list} 
    342  
    343  
    344 % ================================================================ 
    345 % Domain: Horizontal Grid (mesh)  
    346 % ================================================================ 
    347 \section{Horizontal grid mesh (\protect\mdl{domhgr})} 
    348 \label{sec:DOM_hgr} 
    349  
    350 % ------------------------------------------------------------------------------------------------------------- 
    351 %        Coordinates and scale factors  
    352 % ------------------------------------------------------------------------------------------------------------- 
    353 \subsection{Coordinates and scale factors} 
    354 \label{subsec:DOM_hgr_coord_e} 
    355  
    356 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 
    357 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    358 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 
    359 The associated scale factors are defined using the analytical first derivative of the transformation 
    360 \autoref{eq:scale_factors}. 
    361 These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, 
    362 which provide the horizontal and vertical meshes, respectively. 
    363 This section deals with the horizontal mesh parameters. 
    364  
    365 In a horizontal plane, the location of all the model grid points is defined from 
    366 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$. 
    367 The horizontal scale factors are calculated using \autoref{eq:scale_factors}. 
    368 For example, when the longitude and latitude are function of a single value 
    369 ($i$ and $j$, respectively) (geographical configuration of the mesh), 
    370 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 
    371 and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 
    372 The model computes the grid-point positions and scale factors in the horizontal plane as follows: 
     532 
     533Additionally, hybrid combinations of the three main coordinates are available: 
     534$s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e). 
     535 
     536A further choice related to vertical coordinate concerns 
     537the presence (or not) of ocean cavities beneath ice shelves within the model domain. 
     538A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that 
     539the domain contains ocean cavities, 
     540otherwise the top, wet layer of the ocean will always be at the ocean surface. 
     541This option is currently only available for $z$- or $zps$-coordinates. 
     542In the latter case, partial steps are also applied at the ocean/ice shelf interface. 
     543 
     544Within the model, 
     545the arrays describing the grid point depths and vertical scale factors are 
     546three set of three dimensional arrays $(i,j,k)$ defined at 
     547\textit{before}, \textit{now} and \textit{after} time step. 
     548The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
     549They are updated at each model time step. 
     550The initial fixed reference coordinate system is held in variable names with a $\_0$ suffix. 
     551When the linear free surface option is used (\np[=.true.]{ln_linssh}{ln\_linssh}), 
     552\textit{before}, \textit{now} and \textit{after} arrays are initially set to 
     553their reference counterpart and remain fixed. 
     554 
     555%% ================================================================================================= 
     556\subsubsection{Required fields} 
     557\label{sec:DOM_zgr_fields} 
     558 
     559The explicit specification of a range of fields related to the vertical grid are required for 
     560the definition of a configuration. 
     561These include: 
     562 
     563\begin{clines} 
     564int    ln_zco, ln_zps, ln_sco            /* flags for z-coord, z-coord with partial steps and s-coord    */ 
     565int    ln_isfcav                         /* flag  for ice shelf cavities                                 */ 
     566double e3t_1d, e3w_1d                    /* reference vertical scale factors at T and W points           */ 
     567double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ 
     568double e3uw_0, e3vw_0                    /* vertical scale factors 3D coordinate at UW and VW points     */ 
     569int    bottom_level, top_level           /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ 
     570                                         /* For reference:                                               */ 
     571float  bathy_metry                       /* bathymetry used in setting top and bottom levels             */ 
     572\end{clines} 
     573 
     574This set of vertical metrics is sufficient to describe the initial depth and thickness of 
     575every gridcell in the model regardless of the choice of vertical coordinate. 
     576With constant z-levels, e3 metrics will be uniform across each horizontal level. 
     577In the partial step case each e3 at the \jp{bottom\_level} 
     578(and, possibly, \jp{top\_level} if ice cavities are present) 
     579may vary from its horizontal neighbours. 
     580And, in s-coordinates, variations can occur throughout the water column. 
     581With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that 
     582variations occur throughout the water column with displacements related to the sea surface height. 
     583These variations are typically much smaller than those arising from bottom fitted coordinates. 
     584The values for vertical metrics supplied in the domain configuration file can be considered as 
     585those arising from a flat sea surface with zero elevation. 
     586 
     587The \jp{bottom\_level} and \jp{top\_level} 2D arrays define 
     588the \jp{bottom\_level} and top wet levels in each grid column. 
     589Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). 
     590With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. 
     591 
     592%% ================================================================================================= 
     593\subsubsection{Level bathymetry and mask} 
     594\label{subsec:DOM_msk} 
     595 
     596From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: 
    373597\begin{align*} 
    374    \lambda_t &\equiv \text{glamt} =      \lambda (i      ) 
    375   &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\ 
    376    \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2) 
    377   &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\ 
    378    \lambda_v &\equiv \text{glamv} =      \lambda (i      ) 
    379   &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\ 
    380    \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2) 
    381   &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\ 
    382    e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) | 
    383   &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
    384    e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) | 
    385   &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\ 
    386    e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) | 
    387   &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\ 
    388    e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 
    389   &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | 
     598  tmask(i,j,k) &= 
     599  \begin{cases} 
     600    0 &\text{if $                             k <    top\_level(i,j)$} \\ 
     601    1 &\text{if $     bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 
     602    0 &\text{if $k >  bottom\_level(i,j)                            $} 
     603  \end{cases} \\ 
     604  umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     605  vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
     606  fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     607  wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
     608  \text{with~} wmask(i,j,1) &= tmask(i,j,1) 
    390609\end{align*} 
    391 where the last letter of each computational name indicates the grid point considered and 
    392 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 
    393 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 
    394 thus no specific arrays are defined at $w$-points. 
    395  
    396 Note that the definition of the scale factors 
    397 (\ie as the analytical first derivative of the transformation that 
    398 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 
    399 is specific to the \NEMO model \citep{Marti_al_JGR92}. 
    400 As an example, $e_{1t}$ is defined locally at a $t$-point, 
    401 whereas many other models on a C grid choose to define such a scale factor as 
    402 the distance between the $U$-points on each side of the $t$-point. 
    403 Relying on an analytical transformation has two advantages: 
    404 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
    405 since they are first introduced in the continuous equations; 
    406 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
    407 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{Treguier1996}. 
    408 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 
    409 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    410 \begin{figure}[!t] 
    411   \begin{center} 
    412     \includegraphics[]{Fig_zgr_e3} 
    413     \caption{ 
    414       \protect\label{fig:zgr_e3} 
    415       Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
    416       and (b) analytically derived grid-point position and scale factors. 
    417       For both grids here, the same $w$-point depth has been chosen but 
    418       in (a) the $t$-points are set half way between $w$-points while 
    419       in (b) they are defined from an analytical function: 
    420       $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 
    421       Note the resulting difference between the value of the grid-size $\Delta_k$ and 
    422       those of the scale factor $e_k$. 
    423     } 
    424   \end{center} 
    425 \end{figure} 
    426 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    427  
    428 % ------------------------------------------------------------------------------------------------------------- 
    429 %        Choice of horizontal grid 
    430 % ------------------------------------------------------------------------------------------------------------- 
    431 \subsection{Choice of horizontal grid} 
    432 \label{subsec:DOM_hgr_msh_choice} 
    433  
    434 % ------------------------------------------------------------------------------------------------------------- 
    435 %        Grid files 
    436 % ------------------------------------------------------------------------------------------------------------- 
    437 \subsection{Output grid files} 
    438 \label{subsec:DOM_hgr_files} 
    439  
    440 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 
    441 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 
    442 This can be particularly useful for plots and off-line diagnostics. 
    443 In some cases, the user may choose to make a local modification of a scale factor in the code. 
    444 This is the case in global configurations when restricting the width of a specific strait 
    445 (usually a one-grid-point strait that happens to be too wide due to insufficient model resolution). 
    446 An example is Gibraltar Strait in the ORCA2 configuration. 
    447 When such modifications are done, 
    448 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 
    449  
    450 % ================================================================ 
    451 % Domain: Vertical Grid (domzgr) 
    452 % ================================================================ 
    453 \section{Vertical grid (\protect\mdl{domzgr})} 
    454 \label{sec:DOM_zgr} 
    455 %-----------------------------------------nam_zgr & namdom------------------------------------------- 
    456 % 
    457 %\nlst{namzgr}  
    458  
    459 \nlst{namdom}  
    460 %------------------------------------------------------------------------------------------------------------- 
    461  
    462 Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists. 
    463 In the vertical, the model mesh is determined by four things:  
    464 (1) the bathymetry given in meters;  
    465 (2) the number of levels of the model (\jp{jpk});  
    466 (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 
    467 (4) the masking system, \ie the number of wet model levels at each  
    468 $(i,j)$ column of points. 
    469  
    470 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    471 \begin{figure}[!tb] 
    472   \begin{center} 
    473     \includegraphics[]{Fig_z_zps_s_sps} 
    474     \caption{ 
    475       \protect\label{fig:z_zps_s_sps} 
    476       The ocean bottom as seen by the model: 
    477       (a) $z$-coordinate with full step, 
    478       (b) $z$-coordinate with partial step, 
    479       (c) $s$-coordinate: terrain following representation, 
    480       (d) hybrid $s-z$ coordinate, 
    481       (e) hybrid $s-z$ coordinate with partial step, and 
    482       (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). 
    483       Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 
    484     } 
    485   \end{center} 
    486 \end{figure} 
    487 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    488  
    489 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,  
    490 must be done once of all at the beginning of an experiment. 
    491 It is not intended as an option which can be enabled or disabled in the middle of an experiment. 
    492 Three main choices are offered (\autoref{fig:z_zps_s_sps}): 
    493 $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 
    494 $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 
    495 or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 
    496 Hybridation of the three main coordinates are available: 
    497 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 
    498 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 
    499 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 
    500 When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 
    501 the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 
    502 (in other words, the top of the ocean in not a rigid-lid). 
    503 The last choice in terms of vertical coordinate concerns the presence (or not) in 
    504 the model domain of ocean cavities beneath ice shelves. 
    505 Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in. 
    506 This option is currently only available in $z$- or $zps$-coordinate, 
    507 and partial step are also applied at the ocean/ice shelf interface. 
    508  
    509 Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 
    510 reading it from a file. 
    511 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 
    512 \footnote{ 
    513   N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 
    514   so that the computation of the number of wet ocean point in each water column is by-passed}. 
    515 If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 
    516 the ice shelf draft (in meters) is needed. 
    517  
    518 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 
    519 \begin{description} 
    520 \item[\textit{zco}] 
    521   set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 
    522 \item[\textit{zps}] 
    523   set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 
    524   each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 
    525 \item[\textit{sco}] 
    526   smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    527   set the three-dimensional transformation. 
    528 \item[\textit{s-z} and \textit{s-zps}] 
    529   smooth the bathymetry to fulfill the hydrostatic consistency criteria and 
    530   set the three-dimensional transformation $z(i,j,k)$, 
    531   and possibly introduce masking of extra land points to better fit the original bathymetry file. 
    532 \end{description} 
    533 %%% 
    534 \gmcomment{   add the description of the smoothing:  envelop topography...} 
    535 %%% 
    536  
    537 Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 
    538 the arrays describing the grid point depths and vertical scale factors are three set of 
    539 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 
    540 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
    541 They are updated at each model time step using a fixed reference coordinate system which 
    542 computer names have a $\_0$ suffix. 
    543 When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 
    544 \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 
    545  
    546 % ------------------------------------------------------------------------------------------------------------- 
    547 %        Meter Bathymetry 
    548 % ------------------------------------------------------------------------------------------------------------- 
    549 \subsection{Meter bathymetry} 
    550 \label{subsec:DOM_bathy} 
    551  
    552 Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy} 
    553 (found in \ngn{namdom} namelist):  
    554 \begin{description} 
    555 \item[\np{nn\_bathy}~\forcode{= 0}]: 
    556   a flat-bottom domain is defined. 
    557   The total depth $z_w (jpk)$ is given by the coordinate transformation. 
    558   The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 
    559 \item[\np{nn\_bathy}~\forcode{= -1}]: 
    560   a domain with a bump of topography one third of the domain width at the central latitude. 
    561   This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 
    562 \item[\np{nn\_bathy}~\forcode{= 1}]: 
    563   read a bathymetry and ice shelf draft (if needed). 
    564   The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at 
    565   each grid point of the model grid. 
    566   The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto 
    567   the horizontal ocean mesh. 
    568   Defining the bathymetry also defines the coastline: where the bathymetry is zero, 
    569   no model levels are defined (all levels are masked). 
    570  
    571   The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 
    572   each grid point of the model grid. 
    573   This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 
    574   Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 
    575 \end{description} 
    576  
    577 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 
    578 (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 
    579 the rest of the ocean. 
     610 
     611Note that, without ice shelves cavities, 
     612masks at $t-$ and $w-$points are identical with the numerical indexing used 
     613(\autoref{subsec:DOM_Num_Index}). 
     614Nevertheless, 
     615$wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 
     616exactly in the same way as for the bottom boundary. 
     617 
     618%% The specification of closed lateral boundaries requires that at least 
     619%% the first and last rows and columns of the \textit{mbathy} array are set to zero. 
     620%% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 
     621%% the second one and its first column equal to the last but one (and so too the mask arrays) 
     622%% (see \autoref{fig:LBC_jperio}). 
     623 
     624%        Closed seas 
     625%% ================================================================================================= 
     626\subsection{Closed seas} 
     627\label{subsec:DOM_closea} 
     628 
     629When a global ocean is coupled to an atmospheric model it is better to 
     630represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if 
     631the model resolution does not allow their communication with the rest of the ocean. 
    580632This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
    581633so these seas can be removed from the ocean domain. 
    582 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 
    583 but the code has to be adapted to the user's configuration. 
    584  
    585 % ------------------------------------------------------------------------------------------------------------- 
    586 %        z-coordinate  and reference coordinate transformation 
    587 % ------------------------------------------------------------------------------------------------------------- 
    588 \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 
    589             {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 
    590 \label{subsec:DOM_zco} 
    591  
    592 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    593 \begin{figure}[!tb] 
    594   \begin{center} 
    595     \includegraphics[]{Fig_zgr} 
    596     \caption{ 
    597       \protect\label{fig:zgr} 
    598       Default vertical mesh for ORCA2: 30 ocean levels (L30). 
    599       Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from 
    600       \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate. 
    601     } 
    602   \end{center} 
    603 \end{figure} 
    604 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    605  
    606 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 
    607 respectively. 
    608 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 
    609 $gdepw_0(1)$ is the ocean surface. 
    610 There are at most \jp{jpk}-1 $t$-points inside the ocean, 
    611 the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 
    612 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 
    613 analytical derivative with respect to $k$ provides the vertical scale factors. 
    614 The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$. 
    615 This is done in routine \mdl{domzgr} through statement functions, 
    616 using parameters provided in the \ngn{namcfg} namelist. 
    617  
    618 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 
    619 In that case, the parameters \jp{jpk} (number of $w$-levels) and 
    620 \np{pphmax} (total ocean depth in meters) fully define the grid. 
    621  
    622 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 
    623 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps):  
    624 \begin{gather} 
    625   \label{eq:DOM_zgr_ana_1} 
    626     z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 
    627     e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt| 
    628 \end{gather} 
    629 where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 
    630 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 
    631 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 
    632  
    633 If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 
    634 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 
    635 \begin{equation} 
    636   \label{eq:DOM_zgr_ana_2} 
    637   \begin{split} 
    638     e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\ 
    639     e_3^W(k) &= z_T (k    ) - z_T (k - 1) 
    640   \end{split} 
    641 \end{equation} 
    642 This formulation decrease the self-generated circulation into the ice shelf cavity  
    643 (which can, in extreme case, leads to blow up).\\ 
    644   
    645 The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 
    646 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 
    647 This leads to the following conditions: 
    648 \begin{equation} 
    649   \label{eq:DOM_zgr_coef} 
    650   \begin{array}{ll} 
    651     e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\ 
    652     e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 
    653   \end{array} 
    654 \end{equation} 
    655  
    656 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 
    657 the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 
    658 \autoref{eq:DOM_zgr_ana_2} have been determined such that 
    659 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 
    660 For the first standard ORCA2 vertical grid this led to the following values: 
    661 $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 
    662 The resulting depths and scale factors as a function of the model levels are shown in 
    663 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 
    664 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 
    665  
    666 Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 
    667 In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 
    668 in \ngn{namcfg} namelist, and specifies instead the four following parameters: 
    669 \begin{itemize} 
    670 \item 
    671   \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 
    672   The larger \np{ppacr}, the smaller the stretching. 
    673   Values from $3$ to $10$ are usual. 
    674 \item 
    675   \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 
    676   (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 
    677 \item 
    678   \np{ppdzmin}: minimum thickness for the top layer (in meters). 
    679 \item 
    680   \np{pphmax}: total depth of the ocean (meters). 
    681 \end{itemize} 
    682 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 
    683 \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 
    684  
    685 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    686 \begin{table} 
    687   \begin{center} 
    688     \begin{tabular}{c||r|r|r|r} 
    689       \hline 
    690       \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 
    691       \hline 
    692       1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\ 
    693       \hline 
    694       2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\ 
    695       \hline 
    696       3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\ 
    697       \hline 
    698       4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\ 
    699       \hline 
    700       5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\ 
    701       \hline 
    702       6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\ 
    703       \hline 
    704       7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\ 
    705       \hline 
    706       8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\ 
    707       \hline 
    708       9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\ 
    709       \hline 
    710       10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\ 
    711       \hline 
    712       11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\ 
    713       \hline 
    714       12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\ 
    715       \hline 
    716       13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\ 
    717       \hline 
    718       14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\ 
    719       \hline 
    720       15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\ 
    721       \hline 
    722       16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\ 
    723       \hline 
    724       17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\ 
    725       \hline 
    726       18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\ 
    727       \hline 
    728       19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\ 
    729       \hline 
    730       20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\ 
    731       \hline 
    732       21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\ 
    733       \hline 
    734       22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\ 
    735       \hline 
    736       23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\ 
    737       \hline 
    738       24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\ 
    739       \hline 
    740       25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\ 
    741       \hline 
    742       26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\ 
    743       \hline 
    744       27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\ 
    745       \hline 
    746       28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\ 
    747       \hline 
    748       29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\ 
    749       \hline 
    750       30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\ 
    751       \hline 
    752       31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\ 
    753       \hline 
    754     \end{tabular} 
    755   \end{center} 
    756   \caption{ 
    757     \protect\label{tab:orca_zgr} 
    758     Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from 
    759     \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef} 
    760   } 
    761 \end{table} 
    762 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    763  
    764 % ------------------------------------------------------------------------------------------------------------- 
    765 %        z-coordinate with partial step 
    766 % ------------------------------------------------------------------------------------------------------------- 
    767 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 
    768 \label{subsec:DOM_zps} 
    769 %--------------------------------------------namdom------------------------------------------------------- 
    770  
    771 \nlst{namdom}  
    772 %-------------------------------------------------------------------------------------------------------------- 
    773  
    774 In $z$-coordinate partial step, 
    775 the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 
    776 the previous section, \textit{except} in the bottom layer. 
    777 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 
    778 allow a better representation of the bathymetry, especially in the case of small slopes 
    779 (where the bathymetry varies by less than one level thickness from one grid point to the next). 
    780 The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry. 
    781 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 
    782 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 
    783 the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 
    784 This has to be kept in mind when specifying values in \ngn{namdom} namelist, 
    785 as the maximum depth \np{pphmax} in partial steps: 
    786 for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 
    787 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 
    788 Two variables in the namdom namelist are used to define the partial step vertical grid. 
    789 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is 
    790 the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat} 
    791 (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). 
    792  
    793 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  } 
    794  
    795 % ------------------------------------------------------------------------------------------------------------- 
    796 %        s-coordinate 
    797 % ------------------------------------------------------------------------------------------------------------- 
    798 \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 
    799 \label{subsec:DOM_sco} 
    800 %------------------------------------------nam_zgr_sco--------------------------------------------------- 
    801 % 
    802 %\nlst{namzgr_sco}  
    803 %-------------------------------------------------------------------------------------------------------------- 
    804 Options are defined in \ngn{namzgr\_sco}. 
    805 In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 
    806 the product of a depth field and either a stretching function or its derivative, respectively: 
    807  
    808 \begin{align*} 
    809   % \label{eq:DOM_sco_ana} 
    810   z(k)   &= h(i,j) \; z_0 (k) \\ 
    811   e_3(k) &= h(i,j) \; z_0'(k) 
    812 \end{align*} 
    813  
    814 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and 
    815 $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. 
    816 The depth field $h$ is not necessary the ocean depth, 
    817 since a mixed step-like and bottom-following representation of the topography can be used 
    818 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 
    819 The namelist parameter \np{rn\_rmax} determines the slope at which 
    820 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 
    821 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 
    822 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 
    823  
    824 Options for stretching the coordinate are provided as examples, 
    825 but care must be taken to ensure that the vertical stretch used is appropriate for the application. 
    826  
    827 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 
    828 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 
    829 This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: 
    830  
    831 \[ 
    832   z = s_{min} + C (s) (H - s_{min}) 
    833   % \label{eq:SH94_1} 
    834 \] 
    835  
    836 where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 
    837 allows a $z$-coordinate to placed on top of the stretched coordinate, 
    838 and $z$ is the depth (negative down from the asea surface). 
    839 \begin{gather*} 
    840   s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 
    841   % \label{eq:DOM_s} 
    842  \\ 
    843   % \label{eq:DOM_sco_function} 
    844   C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 
    845 \end{gather*} 
    846  
    847 A stretching function, 
    848 modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 
    849 is also available and is more commonly used for shelf seas modelling: 
    850  
    851 \[ 
    852   C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 
    853          + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)} 
    854                         {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)} 
    855   % \label{eq:SH94_2} 
    856 \] 
    857  
    858 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    859 \begin{figure}[!ht] 
    860   \begin{center} 
    861     \includegraphics[]{Fig_sco_function} 
    862     \caption{ 
    863       \protect\label{fig:sco_function} 
    864       Examples of the stretching function applied to a seamount; 
    865       from left to right: surface, surface and bottom, and bottom intensified resolutions 
    866     } 
    867   \end{center} 
    868 \end{figure} 
    869 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    870  
    871 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 
    872 the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 
    873 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 
    874 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution 
    875 (\autoref{fig:sco_function}). 
    876  
    877 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in 
    878 an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. 
    879 In this case the a stretching function $\gamma$ is defined such that: 
    880  
    881 \begin{equation} 
    882   z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 
    883   % \label{eq:z} 
    884 \end{equation} 
    885  
    886 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 
    887  
    888 \begin{gather*} 
    889   % \label{eq:DOM_gamma_deriv} 
    890   \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt) 
    891            + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\ 
    892   \intertext{Where:} 
    893   % \label{eq:DOM_gamma} 
    894   f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 
    895   \quad \text{and} \quad \sigma = \frac{k}{n - 1} 
    896 \end{gather*} 
    897  
    898 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 
    899 the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards 
    900 the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and 
    901 user prescribed surface (\np{rn\_zs}) and bottom depths. 
    902 The bottom cell depth in this example is given as a function of water depth: 
    903  
    904 \[ 
    905   % \label{eq:DOM_zb} 
    906   Z_b = h a + b 
    907 \] 
    908  
    909 where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. 
    910  
    911 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    912 \begin{figure}[!ht] 
    913   \includegraphics[]{Fig_DOM_compare_coordinates_surface} 
    914   \caption{ 
    915     A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), 
    916     a 50 level $Z$-coordinate (contoured surfaces) and 
    917     the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface $100~m$ for 
    918     a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 
    919     For clarity every third coordinate surface is shown. 
    920   } 
    921   \label{fig:fig_compare_coordinates_surface} 
    922 \end{figure} 
    923  % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    924  
    925 This gives a smooth analytical stretching in computational space that is constrained to 
    926 given specified surface and bottom grid cell thicknesses in real space. 
    927 This is not to be confused with the hybrid schemes that 
    928 superimpose geopotential coordinates on terrain following coordinates thus 
    929 creating a non-analytical vertical coordinate that 
    930 therefore may suffer from large gradients in the vertical resolutions. 
    931 This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, 
    932 but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. 
    933  
    934 As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than 
    935 the critical depth $h_c$. 
    936 In this example two options are available in depths shallower than $h_c$, 
    937 with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false 
    938 (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$). 
    939  
    940 Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as 
    941 large slopes lead to hydrostatic consistency. 
    942 A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, 
    943 and is output as part of the model mesh file at the start of the run. 
    944  
    945 % ------------------------------------------------------------------------------------------------------------- 
    946 %        z*- or s*-coordinate 
    947 % ------------------------------------------------------------------------------------------------------------- 
    948 \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 
    949 \label{subsec:DOM_zgr_star} 
    950  
    951 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 
    952  
    953 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances 
    954  
    955 % ------------------------------------------------------------------------------------------------------------- 
    956 %        level bathymetry and mask  
    957 % ------------------------------------------------------------------------------------------------------------- 
    958 \subsection{Level bathymetry and mask} 
    959 \label{subsec:DOM_msk} 
    960  
    961 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 
    962 steps that follow the face of the model cells (step like topography) \citep{Madec_al_JPO96}. 
    963 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 
    964 gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 
    965 mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 
    966 gdept $\leq$ bathy. 
    967  
    968 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after 
    969 mbathy is computed. 
    970 Isolated grid points that do not communicate with another ocean point at the same level are eliminated. 
    971  
    972 As for the representation of bathymetry, a 2D integer array, misfdep, is created. 
    973 misfdep defines the level of the first wet $t$-point. 
    974 All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 
    975 By default, $misfdep(:,:) = 1$ and no cells are masked. 
    976  
    977 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into  
    978 the cavities are performed in the \textit{zgr\_isf} routine. 
    979 The compatibility between ice shelf draft and bathymetry is checked. 
    980 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 
    981 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 
    982 the bathymetry or the ice shelf draft is dug to fit this constrain. 
    983 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 
    984  
    985 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 
    986 \begin{alignat*}{2} 
    987   tmask(i,j,k) &= &  & 
    988     \begin{cases} 
    989                   0 &\text{if $                  k  <    misfdep(i,j)$} \\ 
    990                   1 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\ 
    991                   0 &\text{if $                  k  >     mbathy(i,j)$} 
    992     \end{cases} 
    993   \\ 
    994   umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    995   vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
    996   fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    997                &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    998   wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
    999   \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 
    1000 \end{alignat*} 
    1001  
    1002 Note that, without ice shelves cavities, 
    1003 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 
    1004 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface)  
    1005 exactly in the same way as for the bottom boundary. 
    1006  
    1007 The specification of closed lateral boundaries requires that at least 
    1008 the first and last rows and columns of the \textit{mbathy} array are set to zero. 
    1009 In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 
    1010 the second one and its first column equal to the last but one (and so too the mask arrays) 
    1011 (see \autoref{fig:LBC_jperio}). 
    1012  
    1013 % ================================================================ 
    1014 % Domain: Initial State (dtatsd & istate) 
    1015 % ================================================================ 
    1016 \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    1017 \label{sec:DTA_tsd} 
    1018 %-----------------------------------------namtsd------------------------------------------- 
    1019  
    1020 \nlst{namtsd}  
    1021 %------------------------------------------------------------------------------------------ 
    1022  
    1023 Options are defined in \ngn{namtsd}. 
    1024 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 
    1025 salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 
     634The user has the option to 
     635set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to 
     636optionally decide on the fate of any freshwater imbalance over the area. 
     637The options are explained in \autoref{sec:MISC_closea} but 
     638it should be noted here that a successful use of these options requires 
     639appropriate mask fields to be present in the domain configuration file. 
     640Among the possibilities are: 
     641 
     642\begin{clines} 
     643int closea_mask     /* non-zero values in closed sea areas for optional masking                */ 
     644int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ 
     645int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)   */ 
     646\end{clines} 
     647 
     648%% ================================================================================================= 
     649\subsection{Output grid files} 
     650\label{subsec:DOM_meshmask} 
     651 
     652Most of the arrays relating to a particular ocean model configuration discussed in this chapter 
     653(grid-point position, scale factors) can be saved in a file if 
     654namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to 
     655\forcode{.true.}; 
     656the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. 
     657This is only really useful if 
     658the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and 
     659checking or confirmation is required. 
     660 
     661Alternatively, all the arrays relating to a particular ocean model configuration 
     662(grid-point position, scale factors, depths and masks) can be saved in 
     663a file called \texttt{mesh\_mask} if 
     664namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to 
     665\forcode{.true.}. 
     666This file contains additional fields that can be useful for post-processing applications. 
     667 
     668%% ================================================================================================= 
     669\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
     670\label{sec:DOM_DTA_tsd} 
     671 
     672\begin{listing} 
     673  \nlst{namtsd} 
     674  \caption{\forcode{&namtsd}} 
     675  \label{lst:namtsd} 
     676\end{listing} 
     677 
     678Basic initial state options are defined in \nam{tsd}{tsd}. 
     679By default, the ocean starts from rest (the velocity field is set to zero) and 
     680the initialization of temperature and salinity fields is controlled through the \np{ln_tsd_init}{ln\_tsd\_init} namelist parameter. 
     681 
    1026682\begin{description} 
    1027 \item[\np{ln\_tsd\_init}~\forcode{= .true.}] 
    1028   use a T and S input files that can be given on the model grid itself or on their native input data grid. 
     683\item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on 
     684  the model grid itself or on their native input data grids. 
    1029685  In the latter case, 
    1030686  the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 
    1031687  (see \autoref{subsec:SBC_iof}). 
    1032   The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 
     688  The information relating to the input files are specified in 
     689  the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 
    1033690  The computation is done in the \mdl{dtatsd} module. 
    1034 \item[\np{ln\_tsd\_init}~\forcode{= .false.}] 
    1035   use constant salinity value of $35.5~psu$ and an analytical profile of temperature 
    1036   (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 
     691\item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via 
     692  a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 
     693  The default version sets horizontally uniform T and profiles as used in the GYRE configuration 
     694  (see \autoref{sec:CFGS_gyre}). 
    1037695\end{description} 
    1038696 
    1039 \biblio 
    1040  
    1041 \pindex 
     697\subinc{\input{../../global/epilogue}} 
    1042698 
    1043699\end{document} 

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/doc/latex/NEMO/subfiles/chap_DYN.tex

    r10499 r11830  
    22 
    33\begin{document} 
    4 % ================================================================ 
    5 % Chapter ——— Ocean Dynamics (DYN) 
    6 % ================================================================ 
     4 
    75\chapter{Ocean Dynamics (DYN)} 
    86\label{chap:DYN} 
    97 
    10 \minitoc 
     8\thispagestyle{plain} 
     9 
     10\chaptertoc 
     11 
     12\paragraph{Changes record} ~\\ 
     13 
     14{\footnotesize 
     15  \begin{tabularx}{\textwidth}{l||X|X} 
     16    Release & Author(s) & Modifications \\ 
     17    \hline 
     18    {\em   4.0} & {\em ...} & {\em ...} \\ 
     19    {\em   3.6} & {\em ...} & {\em ...} \\ 
     20    {\em   3.4} & {\em ...} & {\em ...} \\ 
     21    {\em <=3.4} & {\em ...} & {\em ...} 
     22  \end{tabularx} 
     23} 
     24 
     25\clearpage 
    1126 
    1227Using the representation described in \autoref{chap:DOM}, 
    … …  
    3651(surface wind stress calculation using bulk formulae, estimation of mixing coefficients) 
    3752that are carried out in modules SBC, LDF and ZDF and are described in 
    38 \autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.  
     53\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively. 
    3954 
    4055In the present chapter we also describe the diagnostic equations used to compute the horizontal divergence, 
    4156curl of the velocities (\emph{divcur} module) and the vertical velocity (\emph{wzvmod} module). 
    4257 
    43 The different options available to the user are managed by namelist variables.  
    44 For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx},  
     58The different options available to the user are managed by namelist variables. 
     59For term \textit{ttt} in the momentum equations, the logical namelist variables are \textit{ln\_dynttt\_xxx}, 
    4560where \textit{xxx} is a 3 or 4 letter acronym corresponding to each optional scheme. 
    46 If a CPP key is used for this term its name is \key{ttt}. 
     61%If a CPP key is used for this term its name is \key{ttt}. 
    4762The corresponding code can be found in the \textit{dynttt\_xxx} module in the DYN directory, 
    4863and it is usually computed in the \textit{dyn\_ttt\_xxx} subroutine. 
    4964 
    5065The user has the option of extracting and outputting each tendency term from the 3D momentum equations 
    51 (\key{trddyn} defined), as described in \autoref{chap:MISC}. 
    52 Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \key{trdvor} is defined) 
     66(\texttt{trddyn?} defined), as described in \autoref{chap:MISC}. 
     67Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined) 
    5368can be derived from the 3D terms. 
    54 %%% 
    55 \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does  
    56 MISC correspond to "extracting tendency terms" or "vorticity balance"?} 
    57  
    58 % ================================================================ 
    59 % Sea Surface Height evolution & Diagnostics variables 
    60 % ================================================================ 
     69\cmtgm{STEVEN: not quite sure I've got the sense of the last sentence. 
     70  Does MISC correspond to "extracting tendency terms" or "vorticity balance"?} 
     71 
     72%% ================================================================================================= 
    6173\section{Sea surface height and diagnostic variables ($\eta$, $\zeta$, $\chi$, $w$)} 
    6274\label{sec:DYN_divcur_wzv} 
    6375 
    64 %-------------------------------------------------------------------------------------------------------------- 
    65 %           Horizontal divergence and relative vorticity 
    66 %-------------------------------------------------------------------------------------------------------------- 
    67 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
     76%% ================================================================================================= 
     77\subsection[Horizontal divergence and relative vorticity (\textit{divcur.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{divcur})} 
    6878\label{subsec:DYN_divcur} 
    6979 
    70 The vorticity is defined at an $f$-point (\ie corner point) as follows: 
    71 \begin{equation} 
    72   \label{eq:divcur_cur} 
     80The vorticity is defined at an $f$-point (\ie\ corner point) as follows: 
     81\begin{equation} 
     82  \label{eq:DYN_divcur_cur} 
    7383  \zeta =\frac{1}{e_{1f}\,e_{2f} }\left( {\;\delta_{i+1/2} \left[ {e_{2v}\;v} \right] 
    7484      -\delta_{j+1/2} \left[ {e_{1u}\;u} \right]\;} \right) 
    75 \end{equation}  
     85\end{equation} 
    7686 
    7787The horizontal divergence is defined at a $T$-point. 
    7888It is given by: 
    7989\[ 
    80   % \label{eq:divcur_div} 
     90  % \label{eq:DYN_divcur_div} 
    8191  \chi =\frac{1}{e_{1t}\,e_{2t}\,e_{3t} } 
    8292  \left( {\delta_i \left[ {e_{2u}\,e_{3u}\,u} \right] 
    … …  
    96106ensure perfect restartability. 
    97107The vorticity and divergence at the \textit{now} time step are used for the computation of 
    98 the nonlinear advection and of the vertical velocity respectively.  
    99  
    100 %-------------------------------------------------------------------------------------------------------------- 
    101 %           Sea Surface Height evolution 
    102 %-------------------------------------------------------------------------------------------------------------- 
    103 \subsection{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
     108the nonlinear advection and of the vertical velocity respectively. 
     109 
     110%% ================================================================================================= 
     111\subsection[Horizontal divergence and relative vorticity (\textit{sshwzv.F90})]{Horizontal divergence and relative vorticity (\protect\mdl{sshwzv})} 
    104112\label{subsec:DYN_sshwzv} 
    105113 
    106114The sea surface height is given by: 
    107115\begin{equation} 
    108   \label{eq:dynspg_ssh} 
     116  \label{eq:DYN_spg_ssh} 
    109117  \begin{aligned} 
    110118    \frac{\partial \eta }{\partial t} 
    … …  
    115123  \end{aligned} 
    116124\end{equation} 
    117 where \textit{emp} is the surface freshwater budget (evaporation minus precipitation),  
     125where \textit{emp} is the surface freshwater budget (evaporation minus precipitation), 
    118126expressed in Kg/m$^2$/s (which is equal to mm/s), 
    119127and $\rho_w$=1,035~Kg/m$^3$ is the reference density of sea water (Boussinesq approximation). 
    120128If river runoff is expressed as a surface freshwater flux (see \autoref{chap:SBC}) then 
    121 \textit{emp} can be written as the evaporation minus precipitation, minus the river runoff.  
     129\textit{emp} can be written as the evaporation minus precipitation, minus the river runoff. 
    122130The sea-surface height is evaluated using exactly the same time stepping scheme as 
    123 the tracer equation \autoref{eq:tra_nxt}: 
     131the tracer equation \autoref{eq:TRA_nxt}: 
    124132a leapfrog scheme in combination with an Asselin time filter, 
    125 \ie the velocity appearing in \autoref{eq:dynspg_ssh} is centred in time (\textit{now} velocity). 
     133\ie\ the velocity appearing in \autoref{eq:DYN_spg_ssh} is centred in time (\textit{now} velocity). 
    126134This is of paramount importance. 
    127135Replacing $T$ by the number $1$ in the tracer equation and summing over the water column must lead to 
    128136the sea surface height equation otherwise tracer content will not be conserved 
    129 \citep{Griffies_al_MWR01, Leclair_Madec_OM09}. 
     137\citep{griffies.pacanowski.ea_MWR01, leclair.madec_OM09}. 
    130138 
    131139The vertical velocity is computed by an upward integration of the horizontal divergence starting at the bottom, 
    132140taking into account the change of the thickness of the levels: 
    133141\begin{equation} 
    134   \label{eq:wzv} 
     142  \label{eq:DYN_wzv} 
    135143  \left\{ 
    136144    \begin{aligned} 
    … …  
    142150\end{equation} 
    143151 
    144 In the case of a non-linear free surface (\key{vvl}), the top vertical velocity is $-\textit{emp}/\rho_w$,  
     152In the case of a non-linear free surface (\texttt{vvl?}), the top vertical velocity is $-\textit{emp}/\rho_w$, 
    145153as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, 
    146154re-orientated downward. 
    147 \gmcomment{not sure of this...  to be modified with the change in emp setting} 
    148 In the case of a linear free surface, the time derivative in \autoref{eq:wzv} disappears. 
     155\cmtgm{not sure of this...  to be modified with the change in emp setting} 
     156In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. 
    149157The upper boundary condition applies at a fixed level $z=0$. 
    150158The top vertical velocity is thus equal to the divergence of the barotropic transport 
    151 (\ie the first term in the right-hand-side of \autoref{eq:dynspg_ssh}). 
     159(\ie\ the first term in the right-hand-side of \autoref{eq:DYN_spg_ssh}). 
    152160 
    153161Note also that whereas the vertical velocity has the same discrete expression in $z$- and $s$-coordinates, 
    154162its physical meaning is not the same: 
    155163in the second case, $w$ is the velocity normal to the $s$-surfaces. 
    156 Note also that the $k$-axis is re-orientated downwards in the \fortran code compared to 
    157 the indexing used in the semi-discrete equations such as \autoref{eq:wzv} 
    158 (see \autoref{subsec:DOM_Num_Index_vertical}).  
    159  
    160  
    161 % ================================================================ 
    162 % Coriolis and Advection terms: vector invariant form 
    163 % ================================================================ 
     164Note also that the $k$-axis is re-orientated downwards in the \fortran\ code compared to 
     165the indexing used in the semi-discrete equations such as \autoref{eq:DYN_wzv} 
     166(see \autoref{subsec:DOM_Num_Index_vertical}). 
     167 
     168%% ================================================================================================= 
    164169\section{Coriolis and advection: vector invariant form} 
    165170\label{sec:DYN_adv_cor_vect} 
    166 %-----------------------------------------nam_dynadv---------------------------------------------------- 
    167  
    168 \nlst{namdyn_adv}  
    169 %------------------------------------------------------------------------------------------------------------- 
     171 
     172\begin{listing} 
     173  \nlst{namdyn_adv} 
     174  \caption{\forcode{&namdyn_adv}} 
     175  \label{lst:namdyn_adv} 
     176\end{listing} 
    170177 
    171178The vector invariant form of the momentum equations is the one most often used in 
    172 applications of the \NEMO ocean model. 
     179applications of the \NEMO\ ocean model. 
    173180The flux form option (see next section) has been present since version $2$. 
    174 Options are defined through the \ngn{namdyn\_adv} namelist variables Coriolis and 
     181Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables Coriolis and 
    175182momentum advection terms are evaluated using a leapfrog scheme, 
    176 \ie the velocity appearing in these expressions is centred in time (\textit{now} velocity).  
     183\ie\ the velocity appearing in these expressions is centred in time (\textit{now} velocity). 
    177184At the lateral boundaries either free slip, no slip or partial slip boundary conditions are applied following 
    178185\autoref{chap:LBC}. 
    179186 
    180 % ------------------------------------------------------------------------------------------------------------- 
    181 %        Vorticity term  
    182 % ------------------------------------------------------------------------------------------------------------- 
    183 \subsection{Vorticity term (\protect\mdl{dynvor})} 
     187%% ================================================================================================= 
     188\subsection[Vorticity term (\textit{dynvor.F90})]{Vorticity term (\protect\mdl{dynvor})} 
    184189\label{subsec:DYN_vor} 
    185 %------------------------------------------nam_dynvor---------------------------------------------------- 
    186  
    187 \nlst{namdyn_vor}  
    188 %------------------------------------------------------------------------------------------------------------- 
    189  
    190 Options are defined through the \ngn{namdyn\_vor} namelist variables. 
    191 Four discretisations of the vorticity term (\np{ln\_dynvor\_xxx}\forcode{ = .true.}) are available: 
     190 
     191\begin{listing} 
     192  \nlst{namdyn_vor} 
     193  \caption{\forcode{&namdyn_vor}} 
     194  \label{lst:namdyn_vor} 
     195\end{listing} 
     196 
     197Options are defined through the \nam{dyn_vor}{dyn\_vor} namelist variables. 
     198Four discretisations of the vorticity term (\texttt{ln\_dynvor\_xxx}\forcode{=.true.}) are available: 
    192199conserving potential enstrophy of horizontally non-divergent flow (ENS scheme); 
    193200conserving horizontal kinetic energy (ENE scheme); 
    … …  
    195202horizontal kinetic energy for the planetary vorticity term (MIX scheme); 
    196203or conserving both the potential enstrophy of horizontally non-divergent flow and horizontal kinetic energy 
    197 (EEN scheme) (see \autoref{subsec:C_vorEEN}). 
     204(EEN scheme) (see \autoref{subsec:INVARIANTS_vorEEN}). 
    198205In the case of ENS, ENE or MIX schemes the land sea mask may be slightly modified to ensure the consistency of 
    199 vorticity term with analytical equations (\np{ln\_dynvor\_con}\forcode{ = .true.}). 
     206vorticity term with analytical equations (\np[=.true.]{ln_dynvor_con}{ln\_dynvor\_con}). 
    200207The vorticity terms are all computed in dedicated routines that can be found in the \mdl{dynvor} module. 
    201208 
    202 %------------------------------------------------------------- 
    203209%                 enstrophy conserving scheme 
    204 %------------------------------------------------------------- 
    205 \subsubsection{Enstrophy conserving scheme (\protect\np{ln\_dynvor\_ens}\forcode{ = .true.})} 
     210%% ================================================================================================= 
     211\subsubsection[Enstrophy conserving scheme (\forcode{ln_dynvor_ens})]{Enstrophy conserving scheme (\protect\np{ln_dynvor_ens}{ln\_dynvor\_ens})} 
    206212\label{subsec:DYN_vor_ens} 
    207213 
    208214In the enstrophy conserving case (ENS scheme), 
    209215the discrete formulation of the vorticity term provides a global conservation of the enstrophy 
    210 ($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie $\chi$=$0$), 
     216($ [ (\zeta +f ) / e_{3f} ]^2 $ in $s$-coordinates) for a horizontally non-divergent flow (\ie\ $\chi$=$0$), 
    211217but does not conserve the total kinetic energy. 
    212218It is given by: 
    213219\begin{equation} 
    214   \label{eq:dynvor_ens} 
     220  \label{eq:DYN_vor_ens} 
    215221  \left\{ 
    216222    \begin{aligned} 
    … …  
    221227    \end{aligned} 
    222228  \right. 
    223 \end{equation}  
    224  
    225 %------------------------------------------------------------- 
     229\end{equation} 
     230 
    226231%                 energy conserving scheme 
    227 %------------------------------------------------------------- 
    228 \subsubsection{Energy conserving scheme (\protect\np{ln\_dynvor\_ene}\forcode{ = .true.})} 
     232%% ================================================================================================= 
     233\subsubsection[Energy conserving scheme (\forcode{ln_dynvor_ene})]{Energy conserving scheme (\protect\np{ln_dynvor_ene}{ln\_dynvor\_ene})} 
    229234\label{subsec:DYN_vor_ene} 
    230235 
    … …  
    232237It is given by: 
    233238\begin{equation} 
    234   \label{eq:dynvor_ene} 
     239  \label{eq:DYN_vor_ene} 
    235240  \left\{ 
    236241    \begin{aligned} 
    … …  
    241246    \end{aligned} 
    242247  \right. 
    243 \end{equation}  
    244  
    245 %------------------------------------------------------------- 
     248\end{equation} 
     249 
    246250%                 mix energy/enstrophy conserving scheme 
    247 %------------------------------------------------------------- 
    248 \subsubsection{Mixed energy/enstrophy conserving scheme (\protect\np{ln\_dynvor\_mix}\forcode{ = .true.}) } 
     251%% ================================================================================================= 
     252\subsubsection[Mixed energy/enstrophy conserving scheme (\forcode{ln_dynvor_mix})]{Mixed energy/enstrophy conserving scheme (\protect\np{ln_dynvor_mix}{ln\_dynvor\_mix})} 
    249253\label{subsec:DYN_vor_mix} 
    250254 
    251255For the mixed energy/enstrophy conserving scheme (MIX scheme), a mixture of the two previous schemes is used. 
    252 It consists of the ENS scheme (\autoref{eq:dynvor_ens}) for the relative vorticity term, 
    253 and of the ENE scheme (\autoref{eq:dynvor_ene}) applied to the planetary vorticity term. 
    254 \[ 
    255   % \label{eq:dynvor_mix} 
     256It consists of the ENS scheme (\autoref{eq:DYN_vor_ens}) for the relative vorticity term, 
     257and of the ENE scheme (\autoref{eq:DYN_vor_ene}) applied to the planetary vorticity term. 
     258\[ 
     259  % \label{eq:DYN_vor_mix} 
    256260  \left\{ { 
    257261      \begin{aligned} 
    … …  
    268272\] 
    269273 
    270 %------------------------------------------------------------- 
    271274%                 energy and enstrophy conserving scheme 
    272 %------------------------------------------------------------- 
    273 \subsubsection{Energy and enstrophy conserving scheme (\protect\np{ln\_dynvor\_een}\forcode{ = .true.}) } 
     275%% ================================================================================================= 
     276\subsubsection[Energy and enstrophy conserving scheme (\forcode{ln_dynvor_een})]{Energy and enstrophy conserving scheme (\protect\np{ln_dynvor_een}{ln\_dynvor\_een})} 
    274277\label{subsec:DYN_vor_een} 
    275278 
    … …  
    278281the presence of grid point oscillation structures that will be invisible to the operator. 
    279282These structures are \textit{computational modes} that will be at least partly damped by 
    280 the momentum diffusion operator (\ie the subgrid-scale advection), but not by the resolved advection term. 
     283the momentum diffusion operator (\ie\ the subgrid-scale advection), but not by the resolved advection term. 
    281284The ENS and ENE schemes therefore do not contribute to dump any grid point noise in the horizontal velocity field. 
    282285Such noise would result in more noise in the vertical velocity field, an undesirable feature. 
    … …  
    284287$u$ and $v$ are located at different grid points, 
    285288a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. 
    286 \gmcomment{ To circumvent this, Adcroft (ADD REF HERE)  
     289\cmtgm{ To circumvent this, Adcroft (ADD REF HERE) 
    287290Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} 
    288291 
    289 A very nice solution to the problem of double averaging was proposed by \citet{Arakawa_Hsu_MWR90}. 
     292A very nice solution to the problem of double averaging was proposed by \citet{arakawa.hsu_MWR90}. 
    290293The idea is to get rid of the double averaging by considering triad combinations of vorticity. 
    291294It is noteworthy that this solution is conceptually quite similar to the one proposed by 
    292 \citep{Griffies_al_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:C}). 
    293  
    294 The \citet{Arakawa_Hsu_MWR90} vorticity advection scheme for a single layer is modified  
    295 for spherical coordinates as described by \citet{Arakawa_Lamb_MWR81} to obtain the EEN scheme.  
    296 First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point:  
    297 \[ 
    298   % \label{eq:pot_vor} 
     295\citep{griffies.gnanadesikan.ea_JPO98} for the discretization of the iso-neutral diffusion operator (see \autoref{apdx:INVARIANTS}). 
     296 
     297The \citet{arakawa.hsu_MWR90} vorticity advection scheme for a single layer is modified 
     298for spherical coordinates as described by \citet{arakawa.lamb_MWR81} to obtain the EEN scheme. 
     299First consider the discrete expression of the potential vorticity, $q$, defined at an $f$-point: 
     300\[ 
     301  % \label{eq:DYN_pot_vor} 
    299302  q  = \frac{\zeta +f} {e_{3f} } 
    300303\] 
    301 where the relative vorticity is defined by (\autoref{eq:divcur_cur}), 
    302 the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is:  
    303 \begin{equation} 
    304   \label{eq:een_e3f} 
     304where the relative vorticity is defined by (\autoref{eq:DYN_divcur_cur}), 
     305the Coriolis parameter is given by $f=2 \,\Omega \;\sin \varphi _f $ and the layer thickness at $f$-points is: 
     306\begin{equation} 
     307  \label{eq:DYN_een_e3f} 
    305308  e_{3f} = \overline{\overline {e_{3t} }} ^{\,i+1/2,j+1/2} 
    306309\end{equation} 
    307310 
    308 %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    309311\begin{figure}[!ht] 
    310   \begin{center} 
    311     \includegraphics[width=0.70\textwidth]{Fig_DYN_een_triad} 
    312     \caption{ 
    313       \protect\label{fig:DYN_een_triad} 
    314       Triads used in the energy and enstrophy conserving scheme (een) for 
    315       $u$-component (upper panel) and $v$-component (lower panel). 
    316     } 
    317   \end{center} 
     312  \centering 
     313  \includegraphics[width=0.66\textwidth]{DYN_een_triad} 
     314  \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{ 
     315    Triads used in the energy and enstrophy conserving scheme (EEN) for 
     316    $u$-component (upper panel) and $v$-component (lower panel).} 
     317  \label{fig:DYN_een_triad} 
    318318\end{figure} 
    319 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    320  
    321 A key point in \autoref{eq:een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made.  
     319 
     320A key point in \autoref{eq:DYN_een_e3f} is how the averaging in the \textbf{i}- and \textbf{j}- directions is made. 
    322321It uses the sum of masked t-point vertical scale factor divided either by the sum of the four t-point masks 
    323 (\np{nn\_een\_e3f}\forcode{ = 1}), or just by $4$ (\np{nn\_een\_e3f}\forcode{ = .true.}). 
     322(\np[=1]{nn_een_e3f}{nn\_een\_e3f}), or just by $4$ (\np[=.true.]{nn_een_e3f}{nn\_een\_e3f}). 
    324323The latter case preserves the continuity of $e_{3f}$ when one or more of the neighbouring $e_{3t}$ tends to zero and 
    325324extends by continuity the value of $e_{3f}$ into the land areas. 
    … …  
    327326(with a systematic reduction of $e_{3f}$ when a model level intercept the bathymetry) 
    328327that tends to reinforce the topostrophy of the flow 
    329 (\ie the tendency of the flow to follow the isobaths) \citep{Penduff_al_OS07}.  
     328(\ie\ the tendency of the flow to follow the isobaths) \citep{penduff.le-sommer.ea_OS07}. 
    330329 
    331330Next, the vorticity triads, $ {^i_j}\mathbb{Q}^{i_p}_{j_p}$ can be defined at a $T$-point as 
    332331the following triad combinations of the neighbouring potential vorticities defined at f-points 
    333 (\autoref{fig:DYN_een_triad}):  
    334 \begin{equation} 
    335   \label{eq:Q_triads} 
     332(\autoref{fig:DYN_een_triad}): 
     333\begin{equation} 
     334  \label{eq:DYN_Q_triads} 
    336335  _i^j \mathbb{Q}^{i_p}_{j_p} 
    337336  = \frac{1}{12} \ \left(   q^{i-i_p}_{j+j_p} + q^{i+j_p}_{j+i_p} + q^{i+i_p}_{j-j_p}  \right) 
    338337\end{equation} 
    339 where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$.  
    340  
    341 Finally, the vorticity terms are represented as:  
    342 \begin{equation} 
    343   \label{eq:dynvor_een} 
     338where the indices $i_p$ and $k_p$ take the values: $i_p = -1/2$ or $1/2$ and $j_p = -1/2$ or $1/2$. 
     339 
     340Finally, the vorticity terms are represented as: 
     341\begin{equation} 
     342  \label{eq:DYN_vor_een} 
    344343  \left\{ { 
    345344      \begin{aligned} 
    … …  
    350349      \end{aligned} 
    351350    } \right. 
    352 \end{equation}  
     351\end{equation} 
    353352 
    354353This EEN scheme in fact combines the conservation properties of the ENS and ENE schemes. 
    355354It conserves both total energy and potential enstrophy in the limit of horizontally nondivergent flow 
    356 (\ie $\chi$=$0$) (see \autoref{subsec:C_vorEEN}).  
     355(\ie\ $\chi$=$0$) (see \autoref{subsec:INVARIANTS_vorEEN}). 
    357356Applied to a realistic ocean configuration, it has been shown that it leads to a significant reduction of 
    358 the noise in the vertical velocity field \citep{Le_Sommer_al_OM09}. 
     357the noise in the vertical velocity field \citep{le-sommer.penduff.ea_OM09}. 
    359358Furthermore, used in combination with a partial steps representation of bottom topography, 
    360359it improves the interaction between current and topography, 
    361 leading to a larger topostrophy of the flow \citep{Barnier_al_OD06, Penduff_al_OS07}.  
    362  
    363 %-------------------------------------------------------------------------------------------------------------- 
    364 %           Kinetic Energy Gradient term 
    365 %-------------------------------------------------------------------------------------------------------------- 
    366 \subsection{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
     360leading to a larger topostrophy of the flow \citep{barnier.madec.ea_OD06, penduff.le-sommer.ea_OS07}. 
     361 
     362%% ================================================================================================= 
     363\subsection[Kinetic energy gradient term (\textit{dynkeg.F90})]{Kinetic energy gradient term (\protect\mdl{dynkeg})} 
    367364\label{subsec:DYN_keg} 
    368365 
    369 As demonstrated in \autoref{apdx:C}, 
     366As demonstrated in \autoref{apdx:INVARIANTS}, 
    370367there is a single discrete formulation of the kinetic energy gradient term that, 
    371368together with the formulation chosen for the vertical advection (see below), 
    372369conserves the total kinetic energy: 
    373370\[ 
    374   % \label{eq:dynkeg} 
     371  % \label{eq:DYN_keg} 
    375372  \left\{ 
    376373    \begin{aligned} 
    … …  
    381378\] 
    382379 
    383 %-------------------------------------------------------------------------------------------------------------- 
    384 %           Vertical advection term 
    385 %-------------------------------------------------------------------------------------------------------------- 
    386 \subsection{Vertical advection term (\protect\mdl{dynzad}) } 
     380%% ================================================================================================= 
     381\subsection[Vertical advection term (\textit{dynzad.F90})]{Vertical advection term (\protect\mdl{dynzad})} 
    387382\label{subsec:DYN_zad} 
    388383 
    … …  
    391386conserves the total kinetic energy. 
    392387Indeed, the change of KE due to the vertical advection is exactly balanced by 
    393 the change of KE due to the gradient of KE (see \autoref{apdx:C}). 
    394 \[ 
    395   % \label{eq:dynzad} 
     388the change of KE due to the gradient of KE (see \autoref{apdx:INVARIANTS}). 
     389\[ 
     390  % \label{eq:DYN_zad} 
    396391  \left\{ 
    397392    \begin{aligned} 
    … …  
    401396  \right. 
    402397\] 
    403 When \np{ln\_dynzad\_zts}\forcode{ = .true.}, 
     398When \np[=.true.]{ln_dynzad_zts}{ln\_dynzad\_zts}, 
    404399a split-explicit time stepping with 5 sub-timesteps is used on the vertical advection term. 
    405 This option can be useful when the value of the timestep is limited by vertical advection \citep{Lemarie_OM2015}.  
     400This option can be useful when the value of the timestep is limited by vertical advection \citep{lemarie.debreu.ea_OM15}. 
    406401Note that in this case, 
    407402a similar split-explicit time stepping should be used on vertical advection of tracer to ensure a better stability, 
    408 an option which is only available with a TVD scheme (see \np{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 
    409  
    410  
    411 % ================================================================ 
    412 % Coriolis and Advection : flux form 
    413 % ================================================================ 
     403an option which is only available with a TVD scheme (see \np{ln_traadv_tvd_zts}{ln\_traadv\_tvd\_zts} in \autoref{subsec:TRA_adv_tvd}). 
     404 
     405%% ================================================================================================= 
    414406\section{Coriolis and advection: flux form} 
    415407\label{sec:DYN_adv_cor_flux} 
    416 %------------------------------------------nam_dynadv---------------------------------------------------- 
    417  
    418 \nlst{namdyn_adv}  
    419 %------------------------------------------------------------------------------------------------------------- 
    420  
    421 Options are defined through the \ngn{namdyn\_adv} namelist variables. 
     408 
     409Options are defined through the \nam{dyn_adv}{dyn\_adv} namelist variables. 
    422410In the flux form (as in the vector invariant form), 
    423411the Coriolis and momentum advection terms are evaluated using a leapfrog scheme, 
    424 \ie the velocity appearing in their expressions is centred in time (\textit{now} velocity). 
     412\ie\ the velocity appearing in their expressions is centred in time (\textit{now} velocity). 
    425413At the lateral boundaries either free slip, 
    426414no slip or partial slip boundary conditions are applied following \autoref{chap:LBC}. 
    427415 
    428  
    429 %-------------------------------------------------------------------------------------------------------------- 
    430 %           Coriolis plus curvature metric terms 
    431 %-------------------------------------------------------------------------------------------------------------- 
    432 \subsection{Coriolis plus curvature metric terms (\protect\mdl{dynvor}) } 
     416%% ================================================================================================= 
     417\subsection[Coriolis plus curvature metric terms (\textit{dynvor.F90})]{Coriolis plus curvature metric terms (\protect\mdl{dynvor})} 
    433418\label{subsec:DYN_cor_flux} 
    434419 
    435420In flux form, the vorticity term reduces to a Coriolis term in which the Coriolis parameter has been modified to account for the "metric" term. 
    436421This altered Coriolis parameter is thus discretised at $f$-points. 
    437 It is given by:  
     422It is given by: 
    438423\begin{multline*} 
    439   % \label{eq:dyncor_metric} 
     424  % \label{eq:DYN_cor_metric} 
    440425  f+\frac{1}{e_1 e_2 }\left( {v\frac{\partial e_2 }{\partial i}  -  u\frac{\partial e_1 }{\partial j}} \right)  \\ 
    441426  \equiv   f + \frac{1}{e_{1f} e_{2f} } \left( { \ \overline v ^{i+1/2}\delta_{i+1/2} \left[ {e_{2u} } \right] 
    442427      -  \overline u ^{j+1/2}\delta_{j+1/2} \left[ {e_{1u} } \right]  }  \ \right) 
    443 \end{multline*}  
    444  
    445 Any of the (\autoref{eq:dynvor_ens}), (\autoref{eq:dynvor_ene}) and (\autoref{eq:dynvor_een}) schemes can be used to 
     428\end{multline*} 
     429 
     430Any of the (\autoref{eq:DYN_vor_ens}), (\autoref{eq:DYN_vor_ene}) and (\autoref{eq:DYN_vor_een}) schemes can be used to 
    446431compute the product of the Coriolis parameter and the vorticity. 
    447 However, the energy-conserving scheme (\autoref{eq:dynvor_een}) has exclusively been used to date. 
    448 This term is evaluated using a leapfrog scheme, \ie the velocity is centred in time (\textit{now} velocity). 
    449  
    450 %-------------------------------------------------------------------------------------------------------------- 
    451 %           Flux form Advection term 
    452 %-------------------------------------------------------------------------------------------------------------- 
    453 \subsection{Flux form advection term (\protect\mdl{dynadv}) } 
     432However, the energy-conserving scheme (\autoref{eq:DYN_vor_een}) has exclusively been used to date. 
     433This term is evaluated using a leapfrog scheme, \ie\ the velocity is centred in time (\textit{now} velocity). 
     434 
     435%% ================================================================================================= 
     436\subsection[Flux form advection term (\textit{dynadv.F90})]{Flux form advection term (\protect\mdl{dynadv})} 
    454437\label{subsec:DYN_adv_flux} 
    455438 
    456439The discrete expression of the advection term is given by: 
    457440\[ 
    458   % \label{eq:dynadv} 
     441  % \label{eq:DYN_adv} 
    459442  \left\{ 
    460443    \begin{aligned} 
    … …  
    475458a $2^{nd}$ order centered finite difference scheme, CEN2, 
    476459or a $3^{rd}$ order upstream biased scheme, UBS. 
    477 The latter is described in \citet{Shchepetkin_McWilliams_OM05}. 
    478 The schemes are selected using the namelist logicals \np{ln\_dynadv\_cen2} and \np{ln\_dynadv\_ubs}.  
     460The latter is described in \citet{shchepetkin.mcwilliams_OM05}. 
     461The schemes are selected using the namelist logicals \np{ln_dynadv_cen2}{ln\_dynadv\_cen2} and \np{ln_dynadv_ubs}{ln\_dynadv\_ubs}. 
    479462In flux form, the schemes differ by the choice of a space and time interpolation to define the value of 
    480 $u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie at the $T$-, $f$-, 
    481 and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$.  
    482  
    483 %------------------------------------------------------------- 
     463$u$ and $v$ at the centre of each face of $u$- and $v$-cells, \ie\ at the $T$-, $f$-, 
     464and $uw$-points for $u$ and at the $f$-, $T$- and $vw$-points for $v$. 
     465 
    484466%                 2nd order centred scheme 
    485 %------------------------------------------------------------- 
    486 \subsubsection{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln\_dynadv\_cen2}\forcode{ = .true.})} 
     467%% ================================================================================================= 
     468\subsubsection[CEN2: $2^{nd}$ order centred scheme (\forcode{ln_dynadv_cen2})]{CEN2: $2^{nd}$ order centred scheme (\protect\np{ln_dynadv_cen2}{ln\_dynadv\_cen2})} 
    487469\label{subsec:DYN_adv_cen2} 
    488470 
    489471In the centered $2^{nd}$ order formulation, the velocity is evaluated as the mean of the two neighbouring points: 
    490472\begin{equation} 
    491   \label{eq:dynadv_cen2} 
     473  \label{eq:DYN_adv_cen2} 
    492474  \left\{ 
    493475    \begin{aligned} 
    … …  
    496478    \end{aligned} 
    497479  \right. 
    498 \end{equation}  
    499  
    500 The scheme is non diffusive (\ie conserves the kinetic energy) but dispersive (\ie it may create false extrema). 
     480\end{equation} 
     481 
     482The scheme is non diffusive (\ie\ conserves the kinetic energy) but dispersive (\ie\ it may create false extrema). 
    501483It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 
    502484produce a sensible solution. 
    … …  
    504486so $u$ and $v$ are the \emph{now} velocities. 
    505487 
    506 %------------------------------------------------------------- 
    507488%                 UBS scheme 
    508 %------------------------------------------------------------- 
    509 \subsubsection{UBS: Upstream Biased Scheme (\protect\np{ln\_dynadv\_ubs}\forcode{ = .true.})} 
     489%% ================================================================================================= 
     490\subsubsection[UBS: Upstream Biased Scheme (\forcode{ln_dynadv_ubs})]{UBS: Upstream Biased Scheme (\protect\np{ln_dynadv_ubs}{ln\_dynadv\_ubs})} 
    510491\label{subsec:DYN_adv_ubs} 
    511492 
    … …  
    514495For example, the evaluation of $u_T^{ubs} $ is done as follows: 
    515496\begin{equation} 
    516   \label{eq:dynadv_ubs} 
     497  \label{eq:DYN_adv_ubs} 
    517498  u_T^{ubs} =\overline u ^i-\;\frac{1}{6} 
    518499  \begin{cases} 
    … …  
    522503\end{equation} 
    523504where $u"_{i+1/2} =\delta_{i+1/2} \left[ {\delta_i \left[ u \right]} \right]$. 
    524 This results in a dissipatively dominant (\ie hyper-diffusive) truncation error 
    525 \citep{Shchepetkin_McWilliams_OM05}. 
    526 The overall performance of the advection scheme is similar to that reported in \citet{Farrow1995}. 
     505This results in a dissipatively dominant (\ie\ hyper-diffusive) truncation error 
     506\citep{shchepetkin.mcwilliams_OM05}. 
     507The overall performance of the advection scheme is similar to that reported in \citet{farrow.stevens_JPO95}. 
    527508It is a relatively good compromise between accuracy and smoothness. 
    528509It is not a \emph{positive} scheme, meaning that false extrema are permitted. 
    529510But the amplitudes of the false extrema are significantly reduced over those in the centred second order method. 
    530 As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum  
    531 (\ie \np{ln\_dynldf\_lap}\forcode{ = }\np{ln\_dynldf\_bilap}\forcode{ = .false.}), 
     511As the scheme already includes a diffusion component, it can be used without explicit lateral diffusion on momentum 
     512(\ie\ \np[=]{ln_dynldf_lap}{ln\_dynldf\_lap}\np[=.false.]{ln_dynldf_bilap}{ln\_dynldf\_bilap}), 
    532513and it is recommended to do so. 
    533514 
    534515The UBS scheme is not used in all directions. 
    535 In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie $u_{uw}^{ubs}$ and 
    536 $u_{vw}^{ubs}$ in \autoref{eq:dynadv_cen2} are used. 
    537 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the  
     516In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and 
     517$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. 
     518UBS is diffusive and is associated with vertical mixing of momentum. \cmtgm{ gm  pursue the 
    538519sentence:Since vertical mixing of momentum is a source term of the TKE equation...  } 
    539520 
    540 For stability reasons, the first term in (\autoref{eq:dynadv_ubs}), 
     521For stability reasons, the first term in (\autoref{eq:DYN_adv_ubs}), 
    541522which corresponds to a second order centred scheme, is evaluated using the \textit{now} velocity (centred in time), 
    542523while the second term, which is the diffusion part of the scheme, 
    543524is evaluated using the \textit{before} velocity (forward in time). 
    544 This is discussed by \citet{Webb_al_JAOT98} in the context of the Quick advection scheme. 
     525This is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the Quick advection scheme. 
    545526 
    546527Note that the UBS and QUICK (Quadratic Upstream Interpolation for Convective Kinematics) schemes only differ by 
    547528one coefficient. 
    548 Replacing $1/6$ by $1/8$ in (\autoref{eq:dynadv_ubs}) leads to the QUICK advection scheme \citep{Webb_al_JAOT98}. 
     529Replacing $1/6$ by $1/8$ in (\autoref{eq:DYN_adv_ubs}) leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
    549530This option is not available through a namelist parameter, since the $1/6$ coefficient is hard coded. 
    550531Nevertheless it is quite easy to make the substitution in the \mdl{dynadv\_ubs} module and obtain a QUICK scheme. 
    … …  
    553534there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. 
    554535This is an error and should be suppressed soon. 
    555 %%% 
    556 \gmcomment{action :  this have to be done} 
    557 %%% 
    558  
    559 % ================================================================ 
    560 %           Hydrostatic pressure gradient term 
    561 % ================================================================ 
    562 \section{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
     536\cmtgm{action :  this have to be done} 
     537 
     538%% ================================================================================================= 
     539\section[Hydrostatic pressure gradient (\textit{dynhpg.F90})]{Hydrostatic pressure gradient (\protect\mdl{dynhpg})} 
    563540\label{sec:DYN_hpg} 
    564