Ignore:
Timestamp:
2019-11-22T17:15:18+01:00 (18 months ago)
Author:
acc
Message:

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Merge in trunk changes up to 11943 in preparation for end of year merge

Location:
NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc
Files:
15 deleted
65 edited
20 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/README.rst

    r11154 r11954  
     1************************** 
     2Building the documentation 
     3************************** 
    14 
     5.. todo:: 
     6 
     7 
     8 
     9:file:`latex`    : LaTeX sources and Latexmk configuration to build reference manuals with :file:`manual_build.sh` 
     10 
     11:file:`namelists`: Namelist blocks included in the documentation 
     12 
     13:file:`rst`      : |RST man|_ sources and Sphinx configuration to build this guide hereby with :file:`guide_build.sh` 
     14 
     15.. |RST man| replace:: reStructuredText (rst) 
    216 
    317.. warning:: 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/main/bibliography.bib

    r11336 r11954  
    16141614} 
    16151615 
     1616@article{         kraus.turner_tellus67, 
     1617  author = {Kraus, E.B. and Turner, J.}, 
     1618  journal = {Tellus}, 
     1619  pages = {98--106}, 
     1620  title = {A one dimensional model of the seasonal thermocline {II}. {T}he general theory and its consequences}, 
     1621  volume = {19}, 
     1622  year = {1967} 
     1623} 
     1624 
     1625@article{        large.ea_RG97, 
     1626  author        = "Large, W. G. and McWilliams, J. C. and Doney, S. C.", 
     1627  doi           = "10.1029/94RG01872", 
     1628  journal       = "Reviews of Geophysics", 
     1629  number        = {4}, 
     1630  pages         = {363--403}, 
     1631  publisher     = {AGU}, 
     1632  title         = "Oceanic vertical mixing: {A} review and a model with a nonlocal boundary layer parameterization", 
     1633  year          = "1994" 
     1634} 
     1635 
    16161636@techreport{      large.yeager_rpt04, 
    16171637  title         = "Diurnal to decadal global forcing for ocean and sea-ice 
     
    21842204} 
    21852205 
     2206@article{mcwilliams.ea_JFM97, 
     2207   author = {McWilliams, James C. and Sullivan, Peter P. and Moeng, Chin-Hoh}, 
     2208   doi = {10.1017/S0022112096004375}, 
     2209   journal = {Journal of Fluid Mechanics}, 
     2210   pages = {1--30}, 
     2211   title = {Langmuir turbulence in the ocean}, 
     2212   volume = {334}, 
     2213   year = {1997}, 
     2214} 
    21862215@article{         mellor.blumberg_JPO04, 
    21872216  title         = "Wave Breaking and Ocean Surface Layer Thermal Response", 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_DOMAINcfg.tex

    r11599 r11954  
    113113\begin{figure}[!tb] 
    114114  \centering 
    115   \includegraphics[width=0.66\textwidth]{Fig_zgr} 
     115  \includegraphics[width=0.66\textwidth]{DOMCFG_zgr} 
    116116  \caption[DOMAINcfg: default vertical mesh for ORCA2]{ 
    117117    Default vertical mesh for ORCA2: 30 ocean levels (L30). 
     
    444444\begin{figure}[!ht] 
    445445  \centering 
    446   \includegraphics[width=0.66\textwidth]{Fig_sco_function} 
     446  \includegraphics[width=0.66\textwidth]{DOMCFG_sco_function} 
    447447  \caption[DOMAINcfg: examples of the stretching function applied to a seamount]{ 
    448448    Examples of the stretching function applied to a seamount; 
     
    493493\begin{figure}[!ht] 
    494494  \centering 
    495   \includegraphics[width=0.66\textwidth]{Fig_DOM_compare_coordinates_surface} 
     495  \includegraphics[width=0.66\textwidth]{DOMCFG_compare_coordinates_surface} 
    496496  \caption[DOMAINcfg: comparison of $s$- and $z$-coordinate]{ 
    497497    A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines), 
     
    530530This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO\ web site. 
    531531 
    532 \onlyinsubfile{\input{../../global/epilogue}} 
     532\subinc{\input{../../global/epilogue}} 
    533533 
    534534\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_algos.tex

    r11598 r11954  
    311311\begin{figure}[!ht] 
    312312  \centering 
    313   \includegraphics[width=0.66\textwidth]{Fig_ISO_triad} 
     313  %\includegraphics[width=0.66\textwidth]{ALGOS_ISO_triad} 
    314314  \caption[Triads used in the Griffies's like iso-neutral diffision scheme for 
    315315    $u$- and $w$-components)]{ 
     
    461461where $A_{e}$ is the eddy induced velocity coefficient, 
    462462and $r_i$ and $r_j$ the slopes between the iso-neutral and the geopotential surfaces. 
    463 %%gm wrong: to be modified with 2 2D streamfunctions 
     463\cmtgm{Wrong: to be modified with 2 2D streamfunctions} 
    464464In other words, the eddy induced velocity can be derived from a vector streamfuntion, $\phi$, 
    465465which is given by $\phi = A_e\,\textbf{r}$ as $\textbf{U}^*  = \textbf{k} \times \nabla \phi$. 
    466 %%end gm 
    467466 
    468467A traditional way to implement this additional advection is to add it to the eulerian velocity prior to 
     
    822821\ie\ the variance of the tracer is preserved by the discretisation of the skew fluxes. 
    823822 
    824 \onlyinsubfile{\input{../../global/epilogue}} 
     823\subinc{\input{../../global/epilogue}} 
    825824 
    826825\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_diff_opers.tex

    r11598 r11954  
    421421that is a Laplacian diffusion is applied on momentum along the coordinate directions. 
    422422 
    423 \onlyinsubfile{\input{../../global/epilogue}} 
     423\subinc{\input{../../global/epilogue}} 
    424424 
    425425\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_invariants.tex

    r11599 r11954  
    2525\clearpage 
    2626 
    27 %%%  Appendix put in gmcomment as it has not been updated for \zstar and s coordinate 
     27%%%  Appendix put in cmtgm as it has not been updated for \zstar and s coordinate 
    2828%I'm writting this appendix. It will be available in a forthcoming release of the documentation 
    2929 
    30 %\gmcomment{ 
     30%\cmtgm{ 
    3131 
    3232%% ================================================================================================= 
     
    270270 
    271271%gm comment 
    272 \gmcomment{ 
     272\cmtgm{ 
    273273The last equality comes from the following equation, 
    274274\begin{flalign*} 
     
    583583\label{subsec:INVARIANTS_2.6} 
    584584 
    585 \gmcomment{ 
     585\cmtgm{ 
    586586  A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero. 
    587587  In the $z$-coordinate, this property is satisfied locally on a C-grid with 2nd order finite differences 
     
    694694 
    695695%gm comment 
    696 \gmcomment{ 
     696\cmtgm{ 
    697697  \begin{flalign*} 
    698698    \sum\limits_{i,j,k} \biggl\{   p_t\;\partial_t b_t   \biggr\}                                &&&\\ 
     
    14791479%} 
    14801480 
    1481 \onlyinsubfile{\input{../../global/epilogue}} 
     1481\subinc{\input{../../global/epilogue}} 
    14821482 
    14831483\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_s_coord.tex

    r11599 r11954  
    584584the expression of the 3D divergence in the $s-$coordinates established above. 
    585585 
    586 \onlyinsubfile{\input{../../global/epilogue}} 
     586\subinc{\input{../../global/epilogue}} 
    587587 
    588588\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_triads.tex

    r11598 r11954  
    212212\begin{figure}[tb] 
    213213  \centering 
    214   \includegraphics[width=0.66\textwidth]{Fig_GRIFF_triad_fluxes} 
     214  \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_triad_fluxes} 
    215215  \caption[Triads arrangement and tracer gradients to give lateral and vertical tracer fluxes]{ 
    216216    (a) Arrangement of triads $S_i$ and tracer gradients to 
     
    272272\begin{figure}[tb] 
    273273  \centering 
    274   \includegraphics[width=0.66\textwidth]{Fig_GRIFF_qcells} 
     274  \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_qcells} 
    275275  \caption[Triad notation for quarter cells]{ 
    276276    Triad notation for quarter cells. 
     
    657657\begin{figure}[h] 
    658658  \centering 
    659   \includegraphics[width=0.66\textwidth]{Fig_GRIFF_bdry_triads} 
     659  \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_bdry_triads} 
    660660  \caption[Boundary triads]{ 
    661661    (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), 
     
    808808\begin{figure}[h] 
    809809  \centering 
    810   \includegraphics[width=0.66\textwidth]{Fig_GRIFF_MLB_triads} 
     810  \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_MLB_triads} 
    811811  \caption[Definition of mixed-layer depth and calculation of linearly tapered triads]{ 
    812812    Definition of mixed-layer depth and calculation of linearly tapered triads. 
     
    11771177\] 
    11781178 
    1179 \onlyinsubfile{\input{../../global/epilogue}} 
     1179\subinc{\input{../../global/epilogue}} 
    11801180 
    11811181\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_ASM.tex

    r11599 r11954  
    194194\end{clines} 
    195195 
    196 \onlyinsubfile{\input{../../global/epilogue}} 
     196\subinc{\input{../../global/epilogue}} 
    197197 
    198198\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DIA.tex

    r11599 r11954  
    5555A complete description of the use of this I/O server is presented in the next section. 
    5656 
    57 %\gmcomment{                    % start of gmcomment 
     57%\cmtgm{                    % start of gmcomment 
    5858 
    5959%% ================================================================================================= 
     
    19671967\begin{figure}[!t] 
    19681968  \centering 
    1969   \includegraphics[width=0.66\textwidth]{Fig_mask_subasins} 
     1969  \includegraphics[width=0.66\textwidth]{DIA_mask_subasins} 
    19701970  \caption[Decomposition of the World Ocean to compute transports as well as 
    19711971  the meridional stream-function]{ 
     
    20612061The maximum values from the run are also copied to the ocean.output file. 
    20622062 
    2063 \onlyinsubfile{\input{../../global/epilogue}} 
     2063\subinc{\input{../../global/epilogue}} 
    20642064 
    20652065\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DIU.tex

    r11599 r11954  
    160160\] 
    161161 
    162 \onlyinsubfile{\input{../../global/epilogue}} 
     162\subinc{\input{../../global/epilogue}} 
    163163 
    164164\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DOM.tex

    r11598 r11954  
    66\label{chap:DOM} 
    77 
    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.... management of closea sea area : specific to global configuration, both forced and coupled 
    15  
    16 %    {\em 4.0} & {\em Simon M\"{u}ller \& Andrew Coward} & 
    17 %    {\em 
    18 %      Compatibility changes Major simplification has moved many of the options to external domain configuration tools. 
    19 %      (see \autoref{apdx:DOMCFG}) 
    20 %    }                                                                                            \\ 
    21 %    {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} & 
    22 %    {\em First version}                                                                          \\ 
     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 
    2316 
    2417\thispagestyle{plain} 
     
    2922 
    3023{\footnotesize 
    31   \begin{tabularx}{\textwidth}{l||X|X} 
    32     Release & Author(s) & Modifications \\ 
    33     \hline 
    34     {\em   4.0} & {\em ...} & {\em ...} \\ 
    35     {\em   3.6} & {\em ...} & {\em ...} \\ 
    36     {\em   3.4} & {\em ...} & {\em ...} \\ 
    37     {\em <=3.4} & {\em ...} & {\em ...} 
     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                                                                    } 
    3841  \end{tabularx} 
    3942} 
     
    4144\clearpage 
    4245 
    43 Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, 
     46Having defined the continuous equations in \autoref{chap:MB} and 
     47chosen a time discretisation \autoref{chap:TD}, 
    4448we need to choose a grid for spatial discretisation and related numerical algorithms. 
    4549In the present chapter, we provide a general description of the staggered grid used in \NEMO, 
     
    5458\label{subsec:DOM_cell} 
    5559 
    56 \begin{figure}[!tb] 
     60\begin{figure} 
    5761  \centering 
    58   \includegraphics[width=0.66\textwidth]{Fig_cell} 
     62  \includegraphics[width=0.33\textwidth]{DOM_cell} 
    5963  \caption[Arrangement of variables in the unit cell of space domain]{ 
    6064    Arrangement of variables in the unit cell of space domain. 
    6165    $t$ indicates scalar points where 
    6266    temperature, salinity, density, pressure and horizontal divergence are defined. 
    63     $(u,v,w)$ indicates vector points, 
    64     and $f$ indicates vorticity points where 
     67    $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where 
    6568    both relative and planetary vorticities are defined.} 
    6669  \label{fig:DOM_cell} 
    6770\end{figure} 
    6871 
    69 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 
    70 centred second-order finite difference approximation. 
     72The numerical techniques used to solve the Primitive Equations in this model are based on 
     73the traditional, centred second-order finite difference approximation. 
    7174Special attention has been given to the homogeneity of the solution in the three spatial directions. 
    7275The arrangement of variables is the same in all directions. 
    73 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 
    74 the centre of each face of the cells (\autoref{fig:DOM_cell}). 
    75 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 
    76 \citep{mesinger.arakawa_bk76}. 
    77 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 
    78 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 
    79  
    80 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that 
    81 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    82 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. 
    83 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 
    84 the grid-point where the scale factors are defined. 
     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. 
    8591Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 
    8692As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
    8793$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 
    88 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 
    89 while the scale factors are chosen equal to their local analytical value. 
     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. 
    9097An important point here is that the partial derivative of the scale factors must be evaluated by 
    9198centred finite difference approximation, not from their analytical expression. 
    92 This preserves the symmetry of the discrete set of equations and therefore satisfies many of 
    93 the continuous properties (see \autoref{apdx:INVARIANTS}). 
     99This preserves the symmetry of the discrete set of equations and 
     100therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). 
    94101A similar, related remark can be made about the domain size: 
    95 when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors 
    96 (see \autoref{eq:DOM_bar} in the next section). 
    97  
    98 \begin{table}[!tb] 
     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} 
    99106  \centering 
    100   \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 
    101     \hline 
    102     t & $i      $ & $j      $ & $k      $ \\ 
    103     \hline 
    104     u & $i + 1/2$ & $j      $ & $k      $ \\ 
    105     \hline 
    106     v & $i      $ & $j + 1/2$ & $k      $ \\ 
    107     \hline 
    108     w & $i      $ & $j      $ & $k + 1/2$ \\ 
    109     \hline 
    110     f & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
    111     \hline 
    112     uw   & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
    113     \hline 
    114     vw   & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
    115     \hline 
    116     fw   & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
     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$ \\ 
    117124    \hline 
    118125  \end{tabular} 
     
    120127    Location of grid-points as a function of integer or 
    121128    integer and a half value of the column, line or level. 
    122     This indexing is only used for the writing of the semi -discrete equations. 
     129    This indexing is only used for the writing of the semi-discrete equations. 
    123130    In the code, the indexing uses integer values only and 
    124131    is positive downwards in the vertical with $k=1$ at the surface. 
     
    137144firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
    138145since they are first introduced in the continuous equations; 
    139 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
    140 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 
     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}. 
    141150An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 
    142 \begin{figure}[!t] 
     151\begin{figure} 
    143152  \centering 
    144   \includegraphics[width=0.66\textwidth]{Fig_zgr_e3} 
     153  \includegraphics[width=0.5\textwidth]{DOM_zgr_e3} 
    145154  \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ 
    146155    Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
     
    159168\label{subsec:DOM_operators} 
    160169 
    161 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at 
    162 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: 
    163172\begin{alignat*}{2} 
    164173  % \label{eq:DOM_di_mi} 
     
    168177 
    169178Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 
    170 Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has 
    171 its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $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. 
    172183These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
    173 \[ 
     184\begin{gather*} 
    174185  % \label{eq:DOM_grad} 
    175186  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i 
    176187                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j 
    177                   + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 
    178 \] 
    179 \begin{multline*} 
     188                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\ 
    180189  % \label{eq:DOM_lap} 
    181190  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    182191                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) 
    183                             + \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] 
    184193                  + \frac{1}{e_{3t}} 
    185194                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 
    186 \end{multline*} 
    187  
    188 Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
    189 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 
    190200its divergence defined at $t$-points: 
    191 \begin{multline} 
     201\begin{multline*} 
    192202% \label{eq:DOM_curl} 
    193203  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} 
     
    200210                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) 
    201211                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 
    202 \end{multline} 
    203 \begin{equation} 
     212\end{multline*} 
     213\[ 
    204214% \label{eq:DOM_div} 
    205215  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    206216                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 
    207217                              + \frac{1}{e_{3t}} \delta_k (a_3) 
    208 \end{equation} 
    209  
    210 The vertical average over the whole water column is denoted by an overbar and is for 
    211 a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 
     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): 
    212222\begin{equation} 
    213223  \label{eq:DOM_bar} 
     
    215225\end{equation} 
    216226where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 
    217 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over 
    218 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$). 
    219230 
    220231In continuous form, the following properties are satisfied: 
     
    226237\end{gather} 
    227238 
    228 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as 
    229 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 
    230241vector points $(u,v,w)$. 
    231242 
    232243Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 
    233 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
    234 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 
    235 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
    236 \begin{alignat}{4} 
     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} 
    237249  \label{eq:DOM_di_adj} 
    238250  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\ 
     
    241253\end{alignat} 
    242254 
    243 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 
    244257$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 
    245258These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to 
     
    250263\label{subsec:DOM_Num_Index} 
    251264 
    252 \begin{figure}[!tb] 
     265\begin{figure} 
    253266  \centering 
    254   \includegraphics[width=0.66\textwidth]{Fig_index_hor} 
     267  \includegraphics[width=0.33\textwidth]{DOM_index_hor} 
    255268  \caption[Horizontal integer indexing]{ 
    256269    Horizontal integer indexing used in the \fortran\ code. 
     
    261274 
    262275The array representation used in the \fortran\ code requires an integer indexing. 
    263 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 
    264 integer values for $t$-points only while all the other points involve integer and a half values. 
     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. 
    265279Therefore, a specific integer indexing has been defined for points other than $t$-points 
    266280(\ie\ velocity and vorticity grid-points). 
    267 Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. 
     281Furthermore, the direction of the vertical indexing has been reversed and 
     282the surface level set at $k = 1$. 
    268283 
    269284%% ================================================================================================= 
     
    281296\label{subsec:DOM_Num_Index_vertical} 
    282297 
    283 In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in 
    284 the \fortran\ code is the reverse of that used in the semi -discrete equations and 
    285 given in \autoref{subsec:DOM_cell}. 
    286 The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below 
    287 (\autoref{fig:DOM_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}). 
    288303The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 
    289304the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 
    290305Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 
    291306(\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 
    292 in contrast to the indexing on the horizontal plane where the $t$-point has the same index as 
    293 the nearest velocity points in the positive direction of the respective horizontal axis index 
     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 
    294310(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 
    295311Since the scale factors are chosen to be strictly positive, 
     
    298314accommodate the opposing vertical index directions in implementation and documentation. 
    299315 
    300 \begin{figure}[!pt] 
     316\begin{figure} 
    301317  \centering 
    302   \includegraphics[width=0.66\textwidth]{Fig_index_vert} 
     318  \includegraphics[width=0.33\textwidth]{DOM_index_vert} 
    303319  \caption[Vertical integer indexing]{ 
    304320    Vertical integer indexing used in the \fortran\ code. 
     
    314330 
    315331Two typical methods are available to specify the spatial domain configuration; 
    316 they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. 
     332they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in 
     333namelist \nam{cfg}{cfg}. 
    317334 
    318335If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, 
    319 the domain-specific parameters and fields are read from a netCDF input file, 
    320 whose name (without its .nc suffix) can be specified as the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 
     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}. 
    321339 
    322340If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, 
     
    324342subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
    325343These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, 
    326 and default versions that configure the spatial domain for the GYRE reference configuration are present in 
    327 the \path{./src/OCE/USR} directory. 
     344and default versions that configure the spatial domain for the GYRE reference configuration are 
     345present in the \path{./src/OCE/USR} directory. 
    328346 
    329347In version 4.0 there are no longer any options for reading complex bathymetries and 
     
    332350to run similar models with and without partial bottom boxes and/or sigma-coordinates, 
    333351supporting such choices leads to overly complex code. 
    334 Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when 
    335 the model domain itself can be altered by runtime selections. 
    336 The code previously used to perform vertical discretisation has been incorporated into an external tool 
    337 (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 
    338  
    339 The next subsections summarise the parameter and fields related to the configuration of the whole model domain. 
    340 These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code 
    341 inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 
     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. 
    342362The requirements are presented in three sections: 
    343363the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), 
     
    348368\label{subsec:DOM_size} 
    349369 
    350 The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for 
    351 the $i$, $j$ and $k$ directions, respectively. 
    352 Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when 
    353 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 
    354 see \autoref{sec:LBC_mpp}). 
     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}). 
    355375 
    356376The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, 
     
    360380 
    361381The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. 
    362 See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. 
     382See \autoref{sec:LBC_jperio} for details on the available options and 
     383the corresponding values for \jp{jperio}. 
    363384 
    364385%% ================================================================================================= 
     
    370391\label{sec:DOM_hgr_fields} 
    371392 
    372 The explicit specification of a range of mesh-related fields are required for the definition of a configuration. 
     393The explicit specification of a range of mesh-related fields are required for 
     394the definition of a configuration. 
    373395These include: 
    374396 
    375397\begin{clines} 
    376 int    jpiglo, jpjglo, jpkglo            /* global domain sizes                                          */ 
    377 int    jperio                            /* lateral global domain b.c.                                   */ 
    378 double glamt, glamu, glamv, glamf        /* geographic longitude (t,u,v and f points respectively)      */ 
    379 double gphit, gphiu, gphiv, gphif        /* geographic latitude                                          */ 
    380 double e1t, e1u, e1v, e1f                /* horizontal scale factors                                     */ 
    381 double e2t, e2u, e2v, e2f                /* horizontal scale factors                                     */ 
     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                               */ 
    382404\end{clines} 
    383405 
     
    393415 
    394416\begin{clines} 
    395                                          /* Optional:                                                    */ 
    396 int    ORCA, ORCA_index                  /* configuration name, configuration resolution                 */ 
    397 double e1e2u, e1e2v                      /* U and V surfaces (if grid size reduction in some straits)    */ 
    398 double ff_f, ff_t                        /* Coriolis parameter (if not on the sphere)                    */ 
     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)                 */ 
    399421\end{clines} 
    400422 
     
    403425This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points 
    404426(see \autoref{sec:MISC_strait} for illustrated examples). 
    405 The key is to reduce the faces of $T$-cell (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 
     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 
    406429not the volume of the cells. 
    407430Doing otherwise can lead to numerical instability issues. 
    408431In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but 
    409432in cases where a gridsize reduction is required, 
    410 the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or 
    411 pre-computed in \mdl{usrdef\_hgr}. 
    412 If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. 
    413 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set 
    414 the surface-area computation flag: 
     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: 
    415439\texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 
    416440 
    417441\smallskip 
    418442Similar logic applies to the other optional fields: 
    419 \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if 
    420 the mesh is not on a sphere. 
    421 If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 
    422 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set 
    423 the Coriolis computation flag: 
     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: 
    424449\texttt{iff} to a non-zero value to suppress their re-computation. 
    425450 
    426 Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points, 
    427 thus no specific arrays are defined at $w$ points. 
     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. 
    428453 
    429454%% ================================================================================================= 
    430455\subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} 
    431456\label{subsec:DOM_zgr} 
     457 
    432458\begin{listing} 
    433459  \nlst{namdom} 
     
    438464In the vertical, the model mesh is determined by four things: 
    439465\begin{enumerate} 
    440   \item the bathymetry given in meters; 
    441   \item the number of levels of the model (\jp{jpk}); 
    442   \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 
    443   \item the masking system, \ie\ the number of wet model levels at each 
    444 $(i,j)$ location of the horizontal grid. 
     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. 
    445472\end{enumerate} 
    446473 
    447 \begin{figure}[!tb] 
     474\begin{figure} 
    448475  \centering 
    449   \includegraphics[width=0.66\textwidth]{Fig_z_zps_s_sps} 
     476  \includegraphics[width=0.5\textwidth]{DOM_z_zps_s_sps} 
    450477  \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ 
    451478    The ocean bottom as seen by the model: 
    452     (a) $z$-coordinate with full step, 
    453     (b) $z$-coordinate with partial step, 
    454     (c) $s$-coordinate: terrain following representation, 
    455     (d) hybrid $s-z$ coordinate, 
    456     (e) hybrid $s-z$ coordinate with partial step, and 
    457     (f) same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 
    458     Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 
     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).} 
    459489  \label{fig:DOM_z_zps_s_sps} 
    460490\end{figure} 
     
    463493it is not intended to be an option which can be changed in the middle of an experiment. 
    464494The one exception to this statement being the choice of linear or non-linear free surface. 
    465 In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. 
     495In v4.0 the linear free surface option is implemented as 
     496a special case of the non-linear free surface. 
    466497This is computationally wasteful since it uses the structures for time-varying 3D metrics 
    467498for fields that (in the linear free surface case) are fixed. 
    468 However, the linear free-surface is rarely used and implementing it this way means 
    469 a single configuration file can support both options. 
    470  
    471 By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 
    472 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
    473 $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 
    474 When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 
    475 the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface 
     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 
    476510(in other words, the top of the ocean in not a rigid lid). 
    477511 
    478512Note that settings: 
    479 \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} mentioned in the following sections 
    480 appear to be namelist options but they are no longer truly namelist options for \NEMO. 
     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. 
    481516Their value is written to and read from the domain configuration file and 
    482517they should be treated as fixed parameters for a particular configuration. 
    483 They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and 
    484 serve both to provide a record of the choices made whilst building the configuration and 
    485 to trigger appropriate code blocks within \NEMO. 
     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. 
    486521These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. 
    487522 
     
    501536A further choice related to vertical coordinate concerns 
    502537the presence (or not) of ocean cavities beneath ice shelves within the model domain. 
    503 A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, 
     538A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that 
     539the domain contains ocean cavities, 
    504540otherwise the top, wet layer of the ocean will always be at the ocean surface. 
    505541This option is currently only available for $z$- or $zps$-coordinates. 
    506542In the latter case, partial steps are also applied at the ocean/ice shelf interface. 
    507543 
    508 Within the model, the arrays describing the grid point depths and vertical scale factors are three set of 
    509 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 
     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. 
    510548The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
    511549They are updated at each model time step. 
     
    534572\end{clines} 
    535573 
    536 This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in 
    537 the model regardless of the choice of vertical coordinate. 
     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. 
    538576With constant z-levels, e3 metrics will be uniform across each horizontal level. 
    539577In the partial step case each e3 at the \jp{bottom\_level} 
     
    541579may vary from its horizontal neighbours. 
    542580And, in s-coordinates, variations can occur throughout the water column. 
    543 With the non-linear free-surface, all the coordinates behave more like the s-coordinate in 
    544 that variations occur throughout the water column with displacements related to the sea surface height. 
     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. 
    545583These variations are typically much smaller than those arising from bottom fitted coordinates. 
    546584The values for vertical metrics supplied in the domain configuration file can be considered as 
    547585those arising from a flat sea surface with zero elevation. 
    548586 
    549 The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. 
     587The \jp{bottom\_level} and \jp{top\_level} 2D arrays define 
     588the \jp{bottom\_level} and top wet levels in each grid column. 
    550589Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). 
    551590With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. 
     
    556595 
    557596From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: 
    558 \begin{alignat*}{2} 
    559   tmask(i,j,k) &= &  & 
    560     \begin{cases} 
    561                   0 &\text{if $                  k  <    top\_level(i,j)$} \\ 
    562                   1 &\text{if $bottom\_level(i,j) \leq k \leq   top\_level(i,j)$} \\ 
    563                   0 &\text{if $                  k  >     bottom\_level(i,j)$} 
    564     \end{cases} 
    565   \\ 
    566   umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    567   vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
    568   fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    569                &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    570   wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
    571   \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 
    572 \end{alignat*} 
     597\begin{align*} 
     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) 
     609\end{align*} 
    573610 
    574611Note that, without ice shelves cavities, 
    575 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 
    576 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 
     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) 
    577616exactly in the same way as for the bottom boundary. 
    578617 
     
    588627\label{subsec:DOM_closea} 
    589628 
    590 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 
    591 (\eg\ Great Lakes, Caspian sea \dots) even if the model resolution does not allow their communication with 
    592 the rest of the ocean. 
     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. 
    593632This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
    594633so these seas can be removed from the ocean domain. 
    595 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and 
    596 to optionally decide on the fate of any freshwater imbalance over the area. 
    597 The options are explained in \autoref{sec:MISC_closea} but it should be noted here that 
    598 a successful use of these options requires appropriate mask fields to be present in the domain configuration file. 
     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. 
    599640Among the possibilities are: 
    600641 
    601642\begin{clines} 
    602 int    closea_mask          /* non-zero values in closed sea areas for optional masking                  */ 
    603 int    closea_mask_rnf      /* non-zero values in closed sea areas with runoff locations (precip only)  */ 
    604 int    closea_mask_emp      /* non-zero values in closed sea areas with runoff locations (total emp)     */ 
     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)   */ 
    605646\end{clines} 
    606647 
     
    610651 
    611652Most of the arrays relating to a particular ocean model configuration discussed in this chapter 
    612 (grid-point position, scale factors) 
    613 can be saved in a file if 
    614 namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to \forcode{.true.}; 
     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.}; 
    615656the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. 
    616657This is only really useful if 
     
    619660 
    620661Alternatively, all the arrays relating to a particular ocean model configuration 
    621 (grid-point position, scale factors, depths and masks) 
    622 can be saved in a file called \texttt{mesh\_mask} if 
    623 namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. 
     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.}. 
    624666This file contains additional fields that can be useful for post-processing applications. 
    625667 
     
    627669\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    628670\label{sec:DOM_DTA_tsd} 
     671 
    629672\begin{listing} 
    630673  \nlst{namtsd} 
     
    638681 
    639682\begin{description} 
    640 \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on the model grid itself or on their native input data grids. 
    641   In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model 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. 
     685  In the latter case, 
     686  the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 
    642687  (see \autoref{subsec:SBC_iof}). 
    643   The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{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. 
    644690  The computation is done in the \mdl{dtatsd} module. 
    645 \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 
     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}. 
    646693  The default version sets horizontally uniform T and profiles as used in the GYRE configuration 
    647694  (see \autoref{sec:CFGS_gyre}). 
    648695\end{description} 
    649696 
    650 \onlyinsubfile{\input{../../global/epilogue}} 
     697\subinc{\input{../../global/epilogue}} 
    651698 
    652699\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DYN.tex

    r11599 r11954  
    6767Furthermore, the tendency terms associated with the 2D barotropic vorticity balance (when \texttt{trdvor?} is defined) 
    6868can be derived from the 3D terms. 
    69 \gmcomment{STEVEN: not quite sure I've got the sense of the last sentence. does 
    70 MISC correspond to "extracting tendency terms" or "vorticity balance"?} 
     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"?} 
    7171 
    7272%% ================================================================================================= 
     
    153153as changes in the divergence of the barotropic transport are absorbed into the change of the level thicknesses, 
    154154re-orientated downward. 
    155 \gmcomment{not sure of this...  to be modified with the change in emp setting} 
     155\cmtgm{not sure of this...  to be modified with the change in emp setting} 
    156156In the case of a linear free surface, the time derivative in \autoref{eq:DYN_wzv} disappears. 
    157157The upper boundary condition applies at a fixed level $z=0$. 
     
    287287$u$ and $v$ are located at different grid points, 
    288288a price worth paying to avoid a double averaging in the pressure gradient term as in the $B$-grid. 
    289 \gmcomment{ To circumvent this, Adcroft (ADD REF HERE) 
     289\cmtgm{ To circumvent this, Adcroft (ADD REF HERE) 
    290290Nevertheless, this technique strongly distort the phase and group velocity of Rossby waves....} 
    291291 
     
    311311\begin{figure}[!ht] 
    312312  \centering 
    313   \includegraphics[width=0.66\textwidth]{Fig_DYN_een_triad} 
     313  \includegraphics[width=0.66\textwidth]{DYN_een_triad} 
    314314  \caption[Triads used in the energy and enstrophy conserving scheme (EEN)]{ 
    315315    Triads used in the energy and enstrophy conserving scheme (EEN) for 
     
    516516In the vertical, the centred $2^{nd}$ order evaluation of the advection is preferred, \ie\ $u_{uw}^{ubs}$ and 
    517517$u_{vw}^{ubs}$ in \autoref{eq:DYN_adv_cen2} are used. 
    518 UBS is diffusive and is associated with vertical mixing of momentum. \gmcomment{ gm  pursue the 
     518UBS is diffusive and is associated with vertical mixing of momentum. \cmtgm{ gm  pursue the 
    519519sentence:Since vertical mixing of momentum is a source term of the TKE equation...  } 
    520520 
     
    534534there is also the possibility of using a $4^{th}$ order evaluation of the advective velocity as in ROMS. 
    535535This is an error and should be suppressed soon. 
    536 \gmcomment{action :  this have to be done} 
     536\cmtgm{action :  this have to be done} 
    537537 
    538538%% ================================================================================================= 
     
    846846\begin{figure}[!t] 
    847847  \centering 
    848   \includegraphics[width=0.66\textwidth]{Fig_DYN_dynspg_ts} 
     848  \includegraphics[width=0.66\textwidth]{DYN_dynspg_ts} 
    849849  \caption[Split-explicit time stepping scheme for the external and internal modes]{ 
    850850    Schematic of the split-explicit time stepping scheme for the external and internal modes. 
     
    915915it is still significant as shown by \citet{levier.treguier.ea_rpt07} in the case of an analytical barotropic Kelvin wave. 
    916916 
    917 \gmcomment{               %%% copy from griffies Book 
     917\cmtgm{               %%% copy from griffies Book 
    918918 
    919919\textbf{title: Time stepping the barotropic system } 
     
    10431043 
    10441044%% gm %%======>>>>   given here the discrete eqs provided to the solver 
    1045 \gmcomment{               %%% copy from chap-model basics 
     1045\cmtgm{               %%% copy from chap-model basics 
    10461046  \[ 
    10471047    % \label{eq:DYN_spg_flt} 
     
    10541054  and $\mathrm {\mathbf M}$ represents the collected contributions of the Coriolis, hydrostatic pressure gradient, 
    10551055  non-linear and viscous terms in \autoref{eq:MB_dyn}. 
    1056 }   %end gmcomment 
     1056}   %end cmtgm 
    10571057 
    10581058Note that in the linear free surface formulation (\texttt{vvl?} not defined), 
     
    10821082no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}). 
    10831083 
    1084 \gmcomment{ 
     1084\cmtgm{ 
    10851085  Hyperviscous operators are frequently used in the simulation of turbulent flows to 
    10861086  control the dissipation of unresolved small scale features. 
     
    11831183the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen, 
    11841184while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}). 
    1185 \gmcomment{add a remark on the the change in the position of the coefficient} 
     1185\cmtgm{add a remark on the the change in the position of the coefficient} 
    11861186 
    11871187%% ================================================================================================= 
     
    12521252the snow-ice mass is taken into account when computing the surface pressure gradient. 
    12531253 
    1254 \gmcomment{ missing : the lateral boundary condition !!!   another external forcing 
     1254\cmtgm{ missing : the lateral boundary condition !!!   another external forcing 
    12551255 } 
    12561256 
     
    14801480\begin{figure}[!ht] 
    14811481  \centering 
    1482   \includegraphics[width=0.66\textwidth]{Fig_WAD_dynhpg} 
     1482  \includegraphics[width=0.66\textwidth]{DYN_WAD_dynhpg} 
    14831483  \caption[Combinations controlling the limiting of the horizontal pressure gradient in 
    14841484  wetting and drying regimes]{ 
     
    15961596and only array swapping and Asselin filtering is done in \mdl{dynnxt}. 
    15971597 
    1598 \onlyinsubfile{\input{../../global/epilogue}} 
     1598\subinc{\input{../../global/epilogue}} 
    15991599 
    16001600\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_LBC.tex

    r11598 r11954  
    2525\clearpage 
    2626 
    27 %gm% add here introduction to this chapter 
     27\cmtgm{Add here introduction to this chapter} 
    2828 
    2929%% ================================================================================================= 
     
    7979\begin{figure}[!t] 
    8080  \centering 
    81   \includegraphics[width=0.66\textwidth]{Fig_LBC_uv} 
     81  \includegraphics[width=0.66\textwidth]{LBC_uv} 
    8282  \caption[Lateral boundary at $T$-level]{ 
    8383    Lateral boundary (thick line) at T-level. 
     
    104104\begin{figure}[!p] 
    105105  \centering 
    106   \includegraphics[width=0.66\textwidth]{Fig_LBC_shlat} 
     106  \includegraphics[width=0.66\textwidth]{LBC_shlat} 
    107107  \caption[Lateral boundary conditions]{ 
    108108    Lateral boundary conditions 
     
    201201\begin{figure}[!t] 
    202202  \centering 
    203   \includegraphics[width=0.66\textwidth]{Fig_LBC_jperio} 
     203  \includegraphics[width=0.66\textwidth]{LBC_jperio} 
    204204  \caption[Setting of east-west cyclic and symmetric across the Equator boundary conditions]{ 
    205205    Setting of (a) east-west cyclic (b) symmetric across the Equator boundary conditions} 
     
    219219\begin{figure}[!t] 
    220220  \centering 
    221   \includegraphics[width=0.66\textwidth]{Fig_North_Fold_T} 
     221  \includegraphics[width=0.66\textwidth]{LBC_North_Fold_T} 
    222222  \caption[North fold boundary in ORCA 2\deg, 1/4\deg and 1/12\deg]{ 
    223223    North fold boundary with a $T$-point pivot and cyclic east-west boundary condition ($jperio=4$), 
     
    272272\begin{figure}[!t] 
    273273  \centering 
    274   \includegraphics[width=0.66\textwidth]{Fig_mpp} 
     274  \includegraphics[width=0.66\textwidth]{LBC_mpp} 
    275275  \caption{Positioning of a sub-domain when massively parallel processing is used} 
    276276  \label{fig:LBC_mpp} 
     
    325325\begin{figure}[!ht] 
    326326  \centering 
    327   \includegraphics[width=0.66\textwidth]{Fig_mppini2} 
     327  \includegraphics[width=0.66\textwidth]{LBC_mppini2} 
    328328  \caption[Atlantic domain defined for the CLIPPER projet]{ 
    329329    Example of Atlantic domain defined for the CLIPPER projet. 
     
    596596\begin{figure}[!t] 
    597597  \centering 
    598   \includegraphics[width=0.66\textwidth]{Fig_LBC_bdy_geom} 
     598  \includegraphics[width=0.66\textwidth]{LBC_bdy_geom} 
    599599  \caption[Geometry of unstructured open boundary]{Example of geometry of unstructured open boundary} 
    600600  \label{fig:LBC_bdy_geom} 
     
    631631\begin{figure}[!t] 
    632632  \centering 
    633   \includegraphics[width=0.66\textwidth]{Fig_LBC_nc_header} 
     633  \includegraphics[width=0.66\textwidth]{LBC_nc_header} 
    634634  \caption[Header for a \protect\ifile{coordinates.bdy} file]{ 
    635635    Example of the header for a \protect\ifile{coordinates.bdy} file} 
     
    708708direction of rotation). %, e.g. anticlockwise or clockwise. 
    709709 
    710 \onlyinsubfile{\input{../../global/epilogue}} 
     710\subinc{\input{../../global/epilogue}} 
    711711 
    712712\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_LDF.tex

    r11599 r11954  
    6868\label{sec:LDF_slp} 
    6969 
    70 \gmcomment{ 
     70\cmtgm{ 
    7171  we should emphasize here that the implementation is a rather old one. 
    7272  Better work can be achieved by using \citet{griffies.gnanadesikan.ea_JPO98, griffies_bk04} iso-neutral scheme. 
     
    8484$r_{1f}$, $r_{1vw}$, $r_{2t}$, $r_{2vw}$ for $v$. 
    8585 
    86 %gm% add here afigure of the slope in i-direction 
     86\cmtgm{Add here afigure of the slope in i-direction} 
    8787 
    8888%% ================================================================================================= 
     
    9494the diffusive fluxes in the three directions are set to zero and $T$ is assumed to be horizontally uniform, 
    9595\ie\ a linear function of $z_T$, the depth of a $T$-point. 
    96 %gm { Steven : My version is obviously wrong since I'm left with an arbitrary constant which is the local vertical temperature gradient} 
     96\cmtgm{Steven : My version is obviously wrong since 
     97  I'm left with an arbitrary constant which is the local vertical temperature gradient} 
    9798 
    9899\begin{equation} 
     
    112113\end{equation} 
    113114 
    114 %gm%  caution I'm not sure the simplification was a good idea! 
     115\cmtgm{Caution I'm not sure the simplification was a good idea!} 
    115116 
    116117These slopes are computed once in \rou{ldf\_slp\_init} when \np[=.true.]{ln_sco}{ln\_sco}, 
     
    144145\end{equation} 
    145146 
    146 %gm% rewrite this as the explanation is not very clear !!! 
     147\cmtgm{rewrite this as the explanation is not very clear !!!} 
    147148%In practice, \autoref{eq:LDF_slp_iso} is of little help in evaluating the neutral surface slopes. Indeed, for an unsimplified equation of state, the density has a strong dependancy on pressure (here approximated as the depth), therefore applying \autoref{eq:LDF_slp_iso} using the $in situ$ density, $\rho$, computed at T-points leads to a flattening of slopes as the depth increases. This is due to the strong increase of the $in situ$ density with depth. 
    148149 
     
    173174  will include a pressure dependent part, leading to the wrong evaluation of the neutral slopes. 
    174175 
    175 %gm% 
    176176  Note: The solution for $s$-coordinate passes trough the use of different (and better) expression for 
    177177  the constraint on iso-neutral fluxes. 
     
    182182    \alpha \ \textbf{F}(T) = \beta \ \textbf{F}(S) 
    183183  \] 
    184   % gm{  where vector F is ....} 
     184  \cmtgm{where vector F is ....} 
    185185 
    186186This constraint leads to the following definition for the slopes: 
     
    229229This allows an iso-neutral diffusion scheme without additional background horizontal mixing. 
    230230This technique can be viewed as a diffusion operator that acts along large-scale 
    231 (2~$\Delta$x) \gmcomment{2deltax doesnt seem very large scale} iso-neutral surfaces. 
     231(2~$\Delta$x) \cmtgm{2deltax doesnt seem very large scale} iso-neutral surfaces. 
    232232The diapycnal diffusion required for numerical stability is thus minimized and its net effect on the flow is quite small when compared to the effect of an horizontal background mixing. 
    233233 
     
    237237\begin{figure}[!ht] 
    238238  \centering 
    239   \includegraphics[width=0.66\textwidth]{Fig_LDF_ZDF1} 
     239  \includegraphics[width=0.66\textwidth]{LDF_ZDF1} 
    240240  \caption{Averaging procedure for isopycnal slope computation} 
    241241  \label{fig:LDF_ZDF1} 
     
    263263\begin{figure}[!ht] 
    264264  \centering 
    265   \includegraphics[width=0.66\textwidth]{Fig_eiv_slp} 
     265  \includegraphics[width=0.66\textwidth]{LDF_eiv_slp} 
    266266  \caption[Vertical profile of the slope used for lateral mixing in the mixed layer]{ 
    267267    Vertical profile of the slope used for lateral mixing in the mixed layer: 
     
    478478 
    479479%%gm  from Triad appendix  : to be incorporated.... 
    480 \gmcomment{ 
     480\cmtgm{ 
    481481  Values of iso-neutral diffusivity and GM coefficient are set as described in \autoref{sec:LDF_coef}. 
    482482  If none of the keys \key{traldf\_cNd}, N=1,2,3 is set (the default), spatially constant iso-neutral $A_l$ and 
     
    544544\colorbox{yellow}{TBC} 
    545545 
    546 \onlyinsubfile{\input{../../global/epilogue}} 
     546\subinc{\input{../../global/epilogue}} 
    547547 
    548548\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_OBS.tex

    r11598 r11954  
    711711\begin{figure} 
    712712  \centering 
    713   \includegraphics[width=0.66\textwidth]{Fig_OBS_avg_rec} 
     713  \includegraphics[width=0.66\textwidth]{OBS_avg_rec} 
    714714  \caption[Observational weights with a rectangular footprint]{ 
    715715    Weights associated with each model grid box (blue lines and numbers) 
     
    720720\begin{figure} 
    721721  \centering 
    722   \includegraphics[width=0.66\textwidth]{Fig_OBS_avg_rad} 
     722  \includegraphics[width=0.66\textwidth]{OBS_avg_rad} 
    723723  \caption[Observational weights with a radial footprint]{ 
    724724    Weights associated with each model grid box (blue lines and numbers) 
     
    798798\begin{figure} 
    799799  \centering 
    800   \includegraphics[width=0.66\textwidth]{Fig_ASM_obsdist_local} 
     800  \includegraphics[width=0.66\textwidth]{OBS_obsdist_local} 
    801801  \caption[Observations with the geographical distribution]{ 
    802802    Example of the distribution of observations with 
     
    825825\begin{figure} 
    826826  \centering 
    827   \includegraphics[width=0.66\textwidth]{Fig_ASM_obsdist_global} 
     827  \includegraphics[width=0.66\textwidth]{OBS_obsdist_global} 
    828828  \caption[Observations with the round-robin distribution]{ 
    829829    Example of the distribution of observations with 
     
    855855 
    856856%% ================================================================================================= 
    857 \section{Standalone observation operator} 
     857\section{Standalone observation operator (\texttt{SAO})} 
    858858\label{sec:OBS_sao} 
    859859 
     
    11641164\begin{figure} 
    11651165  \centering 
    1166   \includegraphics[width=0.66\textwidth]{Fig_OBS_dataplot_main} 
     1166  \includegraphics[width=0.66\textwidth]{OBS_dataplot_main} 
    11671167  \caption{Main window of dataplot} 
    11681168  \label{fig:OBS_dataplotmain} 
     
    11741174\begin{figure} 
    11751175  \centering 
    1176   \includegraphics[width=0.66\textwidth]{Fig_OBS_dataplot_prof} 
     1176  \includegraphics[width=0.66\textwidth]{OBS_dataplot_prof} 
    11771177  \caption[Profile plot from dataplot]{ 
    11781178    Profile plot from dataplot produced by right clicking on a point in the main window} 
     
    11801180\end{figure} 
    11811181 
    1182 \onlyinsubfile{\input{../../global/epilogue}} 
     1182\subinc{\input{../../global/epilogue}} 
    11831183 
    11841184\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_SBC.tex

    r11599 r11954  
    880880%ENDIF 
    881881 
    882 %\gmcomment{  word doc of runoffs: 
    883 %In the current \NEMO\ setup river runoff is added to emp fluxes, these are then applied at just the sea surface as a volume change (in the variable volume case this is a literal volume change, and in the linear free surface case the free surface is moved) and a salt flux due to the concentration/dilution effect.  There is also an option to increase vertical mixing near river mouths; this gives the effect of having a 3d river.  All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface. 
    884 %Our aim was to code the option to specify the temperature and salinity of river runoff, (as well as the amount), along with the depth that the river water will affect.  This would make it possible to model low salinity outflow, such as the Baltic, and would allow the ocean temperature to be affected by river runoff. 
    885  
    886 %The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between. 
    887  
    888 %To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the tra_sbc module.  We decided to separate them throughout the code, so that the variable emp represented solely evaporation minus precipitation fluxes, and a new 2d variable rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use emp or emps and the changes made are below: 
     882\cmtgm{  word doc of runoffs: 
     883In the current \NEMO\ setup river runoff is added to emp fluxes, 
     884these are then applied at just the sea surface as a volume change (in the variable volume case 
     885this is a literal volume change, and in the linear free surface case the free surface is moved) 
     886and a salt flux due to the concentration/dilution effect. 
     887There is also an option to increase vertical mixing near river mouths; 
     888this gives the effect of having a 3d river. 
     889All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and 
     890at the same temperature as the sea surface. 
     891Our aim was to code the option to specify the temperature and salinity of river runoff, 
     892(as well as the amount), along with the depth that the river water will affect. 
     893This would make it possible to model low salinity outflow, such as the Baltic, 
     894and would allow the ocean temperature to be affected by river runoff. 
     895 
     896The depth option makes it possible to have the river water affecting just the surface layer, 
     897throughout depth, or some specified point in between. 
     898 
     899To do this we need to treat evaporation/precipitation fluxes and river runoff differently in 
     900the \mdl{tra_sbc} module. 
     901We decided to separate them throughout the code, 
     902so that the variable emp represented solely evaporation minus precipitation fluxes, 
     903and a new 2d variable rnf was added which represents the volume flux of river runoff 
     904(in $kg/m^2s$ to remain consistent with $emp$). 
     905This meant many uses of emp and emps needed to be changed, 
     906a list of all modules which use $emp$ or $emps$ and the changes made are below:} 
    889907 
    890908%% ================================================================================================= 
     
    908926  Two different bulk formulae are available: 
    909927 
    910    \begin{description} 
    911    \item [{\np[=1]{nn_isfblk}{nn\_isfblk}}]: The melt rate is based on a balance between the upward ocean heat flux and 
    912      the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}. 
    913    \item [{\np[=2]{nn_isfblk}{nn\_isfblk}}]: The melt rate and the heat flux are based on a 3 equations formulation 
    914      (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation). 
    915      A complete description is available in \citet{jenkins_JGR91}. 
    916    \end{description} 
    917  
    918      Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{losch_JGR08}. 
    919      Its thickness is defined by \np{rn_hisf_tbl}{rn\_hisf\_tbl}. 
    920      The fluxes and friction velocity are computed using the mean temperature, salinity and velocity in the the first \np{rn_hisf_tbl}{rn\_hisf\_tbl} m. 
    921      Then, the fluxes are spread over the same thickness (ie over one or several cells). 
    922      If \np{rn_hisf_tbl}{rn\_hisf\_tbl} larger than top $e_{3}t$, there is no more feedback between the freezing point at the interface and the the top cell temperature. 
    923      This can lead to super-cool temperature in the top cell under melting condition. 
    924      If \np{rn_hisf_tbl}{rn\_hisf\_tbl} smaller than top $e_{3}t$, the top boundary layer thickness is set to the top cell thickness.\\ 
    925  
    926      Each melt bulk formula depends on a exchange coeficient ($\Gamma^{T,S}$) between the ocean and the ice. 
    927      There are 3 different ways to compute the exchange coeficient: 
    928    \begin{description} 
    929         \item [{\np[=0]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are constant and defined by \np{rn_gammas0}{rn\_gammas0} and \np{rn_gammat0}{rn\_gammat0}. 
    930      \begin{gather*} 
     928  \begin{description} 
     929  \item [{\np[=1]{nn_isfblk}{nn\_isfblk}}]: The melt rate is based on a balance between the upward ocean heat flux and 
     930    the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}. 
     931  \item [{\np[=2]{nn_isfblk}{nn\_isfblk}}]: The melt rate and the heat flux are based on a 3 equations formulation 
     932    (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation). 
     933    A complete description is available in \citet{jenkins_JGR91}. 
     934  \end{description} 
     935 
     936  Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{losch_JGR08}. 
     937  Its thickness is defined by \np{rn_hisf_tbl}{rn\_hisf\_tbl}. 
     938  The fluxes and friction velocity are computed using the mean temperature, salinity and velocity in the the first \np{rn_hisf_tbl}{rn\_hisf\_tbl} m. 
     939  Then, the fluxes are spread over the same thickness (ie over one or several cells). 
     940  If \np{rn_hisf_tbl}{rn\_hisf\_tbl} larger than top $e_{3}t$, there is no more feedback between the freezing point at the interface and the the top cell temperature. 
     941  This can lead to super-cool temperature in the top cell under melting condition. 
     942  If \np{rn_hisf_tbl}{rn\_hisf\_tbl} smaller than top $e_{3}t$, the top boundary layer thickness is set to the top cell thickness.\\ 
     943 
     944  Each melt bulk formula depends on a exchange coeficient ($\Gamma^{T,S}$) between the ocean and the ice. 
     945  There are 3 different ways to compute the exchange coeficient: 
     946  \begin{description} 
     947  \item [{\np[=0]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are constant and defined by \np{rn_gammas0}{rn\_gammas0} and \np{rn_gammat0}{rn\_gammat0}. 
     948    \begin{gather*} 
    931949       % \label{eq:SBC_isf_gamma_iso} 
    932        \gamma^{T} = rn\_gammat0 \\ 
    933        \gamma^{S} = rn\_gammas0 
    934      \end{gather*} 
    935      This is the recommended formulation for ISOMIP. 
    936    \item [{\np[=1]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity dependent and defined as 
    937      \begin{gather*} 
    938        \gamma^{T} = rn\_gammat0 \times u_{*} \\ 
    939        \gamma^{S} = rn\_gammas0 \times u_{*} 
    940      \end{gather*} 
    941      where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters). 
    942      See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application. 
    943    \item [{\np[=2]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity and stability dependent and defined as: 
    944 \[ 
    945 \gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}} 
    946 \] 
    947      where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters), 
    948      $\Gamma_{Turb}$ the contribution of the ocean stability and 
    949      $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion. 
    950      See \citet{holland.jenkins_JPO99} for all the details on this formulation. 
    951      This formulation has not been extensively tested in \NEMO\ (not recommended). 
    952    \end{description} 
    953   \item [{\np[=2]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. 
    954    The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting. 
    955    The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL) 
    956    (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front 
    957    (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np[=3]{nn_isf}{nn\_isf}). 
    958    The effective melting length (\np{sn_Leff_isf}{sn\_Leff\_isf}) is read from a file. 
    959   \item [{\np[=3]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. 
    960    The fwf (\np{sn_rnfisf}{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between 
    961    the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and 
    962    the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}). 
    963    The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 
    964   \item [{\np[=4]{nn_isf}{nn\_isf}}]: The ice shelf cavity is opened (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed). 
    965    However, the fwf is not computed but specified from file \np{sn_fwfisf}{sn\_fwfisf}). 
    966    The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 
    967    As in \np[=1]{nn_isf}{nn\_isf}, the fluxes are spread over the top boundary layer thickness (\np{rn_hisf_tbl}{rn\_hisf\_tbl})\\ 
     950      \gamma^{T} = rn\_gammat0 \\ 
     951      \gamma^{S} = rn\_gammas0 
     952    \end{gather*} 
     953    This is the recommended formulation for ISOMIP. 
     954  \item [{\np[=1]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity dependent and defined as 
     955    \begin{gather*} 
     956      \gamma^{T} = rn\_gammat0 \times u_{*} \\ 
     957      \gamma^{S} = rn\_gammas0 \times u_{*} 
     958    \end{gather*} 
     959    where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters). 
     960    See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application. 
     961  \item [{\np[=2]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity and stability dependent and defined as: 
     962    \[ 
     963      \gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}} 
     964    \] 
     965    where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters), 
     966    $\Gamma_{Turb}$ the contribution of the ocean stability and 
     967    $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion. 
     968    See \citet{holland.jenkins_JPO99} for all the details on this formulation. 
     969    This formulation has not been extensively tested in \NEMO\ (not recommended). 
     970  \end{description} 
     971\item [{\np[=2]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. 
     972  The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting. 
     973  The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL) 
     974  (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front 
     975  (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np[=3]{nn_isf}{nn\_isf}). 
     976  The effective melting length (\np{sn_Leff_isf}{sn\_Leff\_isf}) is read from a file. 
     977\item [{\np[=3]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. 
     978  The fwf (\np{sn_rnfisf}{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between 
     979  the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and 
     980  the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}). 
     981  The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 
     982\item [{\np[=4]{nn_isf}{nn\_isf}}]: The ice shelf cavity is opened (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed). 
     983  However, the fwf is not computed but specified from file \np{sn_fwfisf}{sn\_fwfisf}). 
     984  The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 
     985  As in \np[=1]{nn_isf}{nn\_isf}, the fluxes are spread over the top boundary layer thickness (\np{rn_hisf_tbl}{rn\_hisf\_tbl}) 
    968986\end{description} 
    969987 
     
    9861004\begin{figure}[!t] 
    9871005  \centering 
    988   \includegraphics[width=0.66\textwidth]{Fig_SBC_isf} 
     1006  \includegraphics[width=0.66\textwidth]{SBC_isf} 
    9891007  \caption[Ice shelf location and fresh water flux definition]{ 
    9901008    Illustration of the location where the fwf is injected and 
     
    13071325\begin{figure}[!t] 
    13081326  \centering 
    1309   \includegraphics[width=0.66\textwidth]{Fig_SBC_diurnal} 
     1327  \includegraphics[width=0.66\textwidth]{SBC_diurnal} 
    13101328  \caption[Reconstruction of the diurnal cycle variation of short wave flux]{ 
    13111329    Example of reconstruction of the diurnal cycle variation of short wave flux from 
     
    13411359\begin{figure}[!t] 
    13421360  \centering 
    1343   \includegraphics[width=0.66\textwidth]{Fig_SBC_dcy} 
     1361  \includegraphics[width=0.66\textwidth]{SBC_dcy} 
    13441362  \caption[Reconstruction of the diurnal cycle variation of short wave flux on an ORCA2 grid]{ 
    13451363    Example of reconstruction of the diurnal cycle variation of short wave flux from 
     
    15211539% in ocean-ice models. 
    15221540 
    1523 \onlyinsubfile{\input{../../global/epilogue}} 
     1541\subinc{\input{../../global/epilogue}} 
    15241542 
    15251543\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_STO.tex

    r11598 r11954  
    205205The first four parameters define the stochastic part of equation of state. 
    206206 
    207 \onlyinsubfile{\input{../../global/epilogue}} 
     207\subinc{\input{../../global/epilogue}} 
    208208 
    209209\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_TRA.tex

    r11599 r11954  
    1414{\footnotesize 
    1515  \begin{tabularx}{\textwidth}{l||X|X} 
    16     Release & Author(s) & Modifications \\ 
     16    Release          & Author(s)                                   & Modifications      \\ 
    1717    \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 ...} 
     18    {\em        4.0} & {\em Christian \'{E}th\'{e}               } & {\em Review       } \\ 
     19    {\em        3.6} & {\em Gurvan Madec                         } & {\em Update       } \\ 
     20    {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ 
    2221  \end{tabularx} 
    2322} 
     
    3433the tracer equations are available depending on the vertical coordinate used and on the physics used. 
    3534In all the equations presented here, the masking has been omitted for simplicity. 
    36 One must be aware that all the quantities are masked fields and that each time a mean or 
    37 difference operator is used, the resulting field is multiplied by a mask. 
     35One must be aware that all the quantities are masked fields and that 
     36each time a mean or difference operator is used, the resulting field is multiplied by a mask. 
    3837 
    3938The two active tracers are potential temperature and salinity. 
     
    4645NXT stands for next, referring to the time-stepping. 
    4746From left to right, the terms on the rhs of the tracer equations are the advection (ADV), 
    48 the lateral diffusion (LDF), the vertical diffusion (ZDF), the contributions from the external forcings 
    49 (SBC: Surface Boundary Condition, QSR: penetrative Solar Radiation, and BBC: Bottom Boundary Condition), 
    50 the contribution from the bottom boundary Layer (BBL) parametrisation, and an internal damping (DMP) term. 
     47the lateral diffusion (LDF), the vertical diffusion (ZDF), 
     48the contributions from the external forcings (SBC: Surface Boundary Condition, 
     49QSR: penetrative Solar Radiation, and BBC: Bottom Boundary Condition), 
     50the contribution from the bottom boundary Layer (BBL) parametrisation, 
     51and an internal damping (DMP) term. 
    5152The terms QSR, BBC, BBL and DMP are optional. 
    5253The external forcings and parameterisations require complex inputs and complex calculations 
     
    5455LDF and ZDF modules and described in \autoref{chap:SBC}, \autoref{chap:LDF} and 
    5556\autoref{chap:ZDF}, respectively. 
    56 Note that \mdl{tranpc}, the non-penetrative convection module, although located in 
    57 the \path{./src/OCE/TRA} directory as it directly modifies the tracer fields, 
     57Note that \mdl{tranpc}, the non-penetrative convection module, 
     58although located in the \path{./src/OCE/TRA} directory as it directly modifies the tracer fields, 
    5859is described with the model vertical physics (ZDF) together with 
    5960other available parameterization of convection. 
    6061 
    61 In the present chapter we also describe the diagnostic equations used to compute the sea-water properties 
    62 (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and freezing point with 
    63 associated modules \mdl{eosbn2} and \mdl{phycst}). 
     62In the present chapter we also describe the diagnostic equations used to 
     63compute the sea-water properties (density, Brunt-V\"{a}is\"{a}l\"{a} frequency, specific heat and 
     64freezing point with associated modules \mdl{eosbn2} and \mdl{phycst}). 
    6465 
    6566The different options available to the user are managed by namelist logicals. 
     
    7071 
    7172The user has the option of extracting each tendency term on the RHS of the tracer equation for output 
    72 (\np{ln_tra_trd}{ln\_tra\_trd} or \np[=.true.]{ln_tra_mxl}{ln\_tra\_mxl}), as described in \autoref{chap:DIA}. 
     73(\np{ln_tra_trd}{ln\_tra\_trd} or \np[=.true.]{ln_tra_mxl}{ln\_tra\_mxl}), 
     74as described in \autoref{chap:DIA}. 
    7375 
    7476%% ================================================================================================= 
     
    8587the advection tendency of a tracer is expressed in flux form, 
    8688\ie\ as the divergence of the advective fluxes. 
    87 Its discrete expression is given by : 
     89Its discrete expression is given by: 
    8890\begin{equation} 
    8991  \label{eq:TRA_adv} 
     
    9496where $\tau$ is either T or S, and $b_t = e_{1t} \, e_{2t} \, e_{3t}$ is the volume of $T$-cells. 
    9597The flux form in \autoref{eq:TRA_adv} implicitly requires the use of the continuity equation. 
    96 Indeed, it is obtained by using the following equality: $\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which 
    97 results from the use of the continuity equation, $\partial_t e_3 + e_3 \; \nabla \cdot \vect U = 0$ 
    98 (which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, \ie\ \np[=.true.]{ln_linssh}{ln\_linssh}). 
    99 Therefore it is of paramount importance to design the discrete analogue of the advection tendency so that 
    100 it is consistent with the continuity equation in order to enforce the conservation properties of 
    101 the continuous equations. 
    102 In other words, by setting $\tau = 1$ in (\autoref{eq:TRA_adv}) we recover the discrete form of 
    103 the continuity equation which is used to calculate the vertical velocity. 
    104 \begin{figure}[!t] 
     98Indeed, it is obtained by using the following equality: 
     99$\nabla \cdot (\vect U \, T) = \vect U \cdot \nabla T$ which 
     100results from the use of the continuity equation, 
     101$\partial_t e_3 + e_3 \; \nabla \cdot \vect U = 0$ 
     102(which reduces to $\nabla \cdot \vect U = 0$ in linear free surface, 
     103\ie\ \np[=.true.]{ln_linssh}{ln\_linssh}). 
     104Therefore it is of paramount importance to 
     105design the discrete analogue of the advection tendency so that 
     106it is consistent with the continuity equation in order to 
     107enforce the conservation properties of the continuous equations. 
     108In other words, by setting $\tau = 1$ in (\autoref{eq:TRA_adv}) we recover 
     109the discrete form of the continuity equation which is used to calculate the vertical velocity. 
     110\begin{figure} 
    105111  \centering 
    106   \includegraphics[width=0.66\textwidth]{Fig_adv_scheme} 
     112  \includegraphics[width=0.66\textwidth]{TRA_adv_scheme} 
    107113  \caption[Ways to evaluate the tracer value and the amount of tracer exchanged]{ 
    108114    Schematic representation of some ways used to evaluate the tracer value at $u$-point and 
     
    120126\end{figure} 
    121127 
    122 The key difference between the advection schemes available in \NEMO\ is the choice made in space and 
    123 time interpolation to define the value of the tracer at the velocity points 
     128The key difference between the advection schemes available in \NEMO\ is the choice made in 
     129space and time interpolation to define the value of the tracer at the velocity points 
    124130(\autoref{fig:TRA_adv_scheme}). 
    125131 
     
    129135 
    130136\begin{description} 
    131 \item [linear free surface:] (\np[=.true.]{ln_linssh}{ln\_linssh}) 
     137\item [linear free surface] (\np[=.true.]{ln_linssh}{ln\_linssh}) 
    132138  the first level thickness is constant in time: 
    133   the vertical boundary condition is applied at the fixed surface $z = 0$ rather than on 
    134   the moving surface $z = \eta$. 
    135   There is a non-zero advective flux which is set for all advection schemes as 
    136   $\tau_w|_{k = 1/2} = T_{k = 1}$, \ie\ the product of surface velocity (at $z = 0$) by 
    137   the first level tracer value. 
    138 \item [non-linear free surface:] (\np[=.false.]{ln_linssh}{ln\_linssh}) 
     139  the vertical boundary condition is applied at the fixed surface $z = 0$ rather than 
     140  on the moving surface $z = \eta$. 
     141  There is a non-zero advective flux which is set for 
     142  all advection schemes as $\tau_w|_{k = 1/2} = T_{k = 1}$, 
     143  \ie\ the product of surface velocity (at $z = 0$) by the first level tracer value. 
     144\item [non-linear free surface] (\np[=.false.]{ln_linssh}{ln\_linssh}) 
    139145  convergence/divergence in the first ocean level moves the free surface up/down. 
    140   There is no tracer advection through it so that the advective fluxes through the surface are also zero. 
     146  There is no tracer advection through it so that 
     147  the advective fluxes through the surface are also zero. 
    141148\end{description} 
    142149 
    143150In all cases, this boundary condition retains local conservation of tracer. 
    144 Global conservation is obtained in non-linear free surface case, but \textit{not} in the linear free surface case. 
    145 Nevertheless, in the latter case, it is achieved to a good approximation since 
    146 the non-conservative term is the product of the time derivative of the tracer and the free surface height, 
    147 two quantities that are not correlated \citep{roullet.madec_JGR00, griffies.pacanowski.ea_MWR01, campin.adcroft.ea_OM04}. 
     151Global conservation is obtained in non-linear free surface case, 
     152but \textit{not} in the linear free surface case. 
     153Nevertheless, in the latter case, 
     154it is achieved to a good approximation since the non-conservative term is 
     155the product of the time derivative of the tracer and the free surface height, 
     156two quantities that are not correlated 
     157\citep{roullet.madec_JGR00, griffies.pacanowski.ea_MWR01, campin.adcroft.ea_OM04}. 
    148158 
    149159The velocity field that appears in (\autoref{eq:TRA_adv} is 
     
    153163(see \autoref{chap:LDF}). 
    154164 
    155 Several tracer advection scheme are proposed, namely a $2^{nd}$ or $4^{th}$ order centred schemes (CEN), 
    156 a $2^{nd}$ or $4^{th}$ order Flux Corrected Transport scheme (FCT), a Monotone Upstream Scheme for 
    157 Conservative Laws scheme (MUSCL), a $3^{rd}$ Upstream Biased Scheme (UBS, also often called UP3), 
    158 and a Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms scheme (QUICKEST). 
    159 The choice is made in the \nam{tra_adv}{tra\_adv} namelist, by setting to \forcode{.true.} one of 
    160 the logicals \textit{ln\_traadv\_xxx}. 
    161 The corresponding code can be found in the \textit{traadv\_xxx.F90} module, where 
    162 \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 
    163 By default (\ie\ in the reference namelist, \textit{namelist\_ref}), all the logicals are set to \forcode{.false.}. 
    164 If the user does not select an advection scheme in the configuration namelist (\textit{namelist\_cfg}), 
    165 the tracers will \textit{not} be advected! 
     165Several tracer advection scheme are proposed, 
     166namely a $2^{nd}$ or $4^{th}$ order \textbf{CEN}tred schemes (CEN), 
     167a $2^{nd}$ or $4^{th}$ order \textbf{F}lux \textbf{C}orrected \textbf{T}ransport scheme (FCT), 
     168a \textbf{M}onotone \textbf{U}pstream \textbf{S}cheme for 
     169\textbf{C}onservative \textbf{L}aws scheme (MUSCL), 
     170a $3^{rd}$ \textbf{U}pstream \textbf{B}iased \textbf{S}cheme (UBS, also often called UP3), 
     171and a \textbf{Q}uadratic \textbf{U}pstream \textbf{I}nterpolation for 
     172\textbf{C}onvective \textbf{K}inematics with 
     173\textbf{E}stimated \textbf{S}treaming \textbf{T}erms scheme (QUICKEST). 
     174The choice is made in the \nam{tra_adv}{tra\_adv} namelist, 
     175by setting to \forcode{.true.} one of the logicals \textit{ln\_traadv\_xxx}. 
     176The corresponding code can be found in the \textit{traadv\_xxx.F90} module, 
     177where \textit{xxx} is a 3 or 4 letter acronym corresponding to each scheme. 
     178By default (\ie\ in the reference namelist, \textit{namelist\_ref}), 
     179all the logicals are set to \forcode{.false.}. 
     180If the user does not select an advection scheme in the configuration namelist 
     181(\textit{namelist\_cfg}), the tracers will \textit{not} be advected! 
    166182 
    167183Details of the advection schemes are given below. 
    168 The choosing an advection scheme is a complex matter which depends on the model physics, model resolution, 
    169 type of tracer, as well as the issue of numerical cost. In particular, we note that 
     184The choosing an advection scheme is a complex matter which depends on the 
     185model physics, model resolution, type of tracer, as well as the issue of numerical cost. 
     186In particular, we note that 
    170187 
    171188\begin{enumerate} 
    172 \item CEN and FCT schemes require an explicit diffusion operator while the other schemes are diffusive enough so that 
    173   they do not necessarily need additional diffusion; 
    174 \item CEN and UBS are not \textit{positive} schemes 
    175   \footnote{negative values can appear in an initially strictly positive tracer field which is advected}, 
     189\item CEN and FCT schemes require an explicit diffusion operator while 
     190  the other schemes are diffusive enough so that they do not necessarily need additional diffusion; 
     191\item CEN and UBS are not \textit{positive} schemes \footnote{negative values can appear in 
     192    an initially strictly positive tracer field which is advected}, 
    176193  implying that false extrema are permitted. 
    177194  Their use is not recommended on passive tracers; 
    178 \item It is recommended that the same advection-diffusion scheme is used on both active and passive tracers. 
     195\item It is recommended that the same advection-diffusion scheme is used on 
     196  both active and passive tracers. 
    179197\end{enumerate} 
    180198 
    181 Indeed, if a source or sink of a passive tracer depends on an active one, the difference of treatment of active and 
    182 passive tracers can create very nice-looking frontal structures that are pure numerical artefacts. 
     199Indeed, if a source or sink of a passive tracer depends on an active one, 
     200the difference of treatment of active and passive tracers can create 
     201very nice-looking frontal structures that are pure numerical artefacts. 
    183202Nevertheless, most of our users set a different treatment on passive and active tracers, 
    184203that's the reason why this possibility is offered. 
    185 We strongly suggest them to perform a sensitivity experiment using a same treatment to assess the robustness of 
    186 their results. 
     204We strongly suggest them to perform a sensitivity experiment using a same treatment to 
     205assess the robustness of their results. 
    187206 
    188207%% ================================================================================================= 
     
    192211%        2nd order centred scheme 
    193212 
    194 The centred advection scheme (CEN) is used when \np[=.true.]{ln_traadv_cen}{ln\_traadv\_cen}. 
    195 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 
     213The \textbf{CEN}tred advection scheme (CEN) is used when \np[=.true.]{ln_traadv_cen}{ln\_traadv\_cen}. 
     214Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on 
     215horizontal (iso-level) and vertical direction by 
    196216setting \np{nn_cen_h}{nn\_cen\_h} and \np{nn_cen_v}{nn\_cen\_v} to $2$ or $4$. 
    197217CEN implementation can be found in the \mdl{traadv\_cen} module. 
    198218 
    199 In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points is evaluated as the mean of 
    200 the two neighbouring $T$-point values. 
     219In the $2^{nd}$ order centred formulation (CEN2), the tracer at velocity points is evaluated as 
     220the mean of the two neighbouring $T$-point values. 
    201221For example, in the $i$-direction : 
    202222\begin{equation} 
     
    205225\end{equation} 
    206226 
    207 CEN2 is non diffusive (\ie\ it conserves the tracer variance, $\tau^2$) but dispersive 
    208 (\ie\ it may create false extrema). 
    209 It is therefore notoriously noisy and must be used in conjunction with an explicit diffusion operator to 
    210 produce a sensible solution. 
    211 The associated time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 
     227CEN2 is non diffusive (\ie\ it conserves the tracer variance, $\tau^2$) but 
     228dispersive (\ie\ it may create false extrema). 
     229It is therefore notoriously noisy and must be used in conjunction with 
     230an explicit diffusion operator to produce a sensible solution. 
     231The associated time-stepping is performed using 
     232a leapfrog scheme in conjunction with an Asselin time-filter, 
    212233so $T$ in (\autoref{eq:TRA_adv_cen2}) is the \textit{now} tracer value. 
    213234 
     
    217238%        4nd order centred scheme 
    218239 
    219 In the $4^{th}$ order formulation (CEN4), tracer values are evaluated at u- and v-points as 
    220 a $4^{th}$ order interpolation, and thus depend on the four neighbouring $T$-points. 
     240In the $4^{th}$ order formulation (CEN4), 
     241tracer values are evaluated at u- and v-points as a $4^{th}$ order interpolation, 
     242and thus depend on the four neighbouring $T$-points. 
    221243For example, in the $i$-direction: 
    222244\begin{equation} 
     
    226248In the vertical direction (\np[=4]{nn_cen_v}{nn\_cen\_v}), 
    227249a $4^{th}$ COMPACT interpolation has been prefered \citep{demange_phd14}. 
    228 In the COMPACT scheme, both the field and its derivative are interpolated, which leads, after a matrix inversion, 
    229 spectral characteristics similar to schemes of higher order \citep{lele_JCP92}. 
     250In the COMPACT scheme, both the field and its derivative are interpolated, 
     251which leads, after a matrix inversion, spectral characteristics similar to schemes of higher order 
     252\citep{lele_JCP92}. 
    230253 
    231254Strictly speaking, the CEN4 scheme is not a $4^{th}$ order advection scheme but 
    232255a $4^{th}$ order evaluation of advective fluxes, 
    233256since the divergence of advective fluxes \autoref{eq:TRA_adv} is kept at $2^{nd}$ order. 
    234 The expression \textit{$4^{th}$ order scheme} used in oceanographic literature is usually associated with 
    235 the scheme presented here. 
    236 Introducing a \forcode{.true.} $4^{th}$ order advection scheme is feasible but, for consistency reasons, 
    237 it requires changes in the discretisation of the tracer advection together with changes in the continuity equation, 
    238 and the momentum advection and pressure terms. 
     257The expression \textit{$4^{th}$ order scheme} used in oceanographic literature is 
     258usually associated with the scheme presented here. 
     259Introducing a ``true'' $4^{th}$ order advection scheme is feasible but, for consistency reasons, 
     260it requires changes in the discretisation of the tracer advection together with 
     261changes in the continuity equation, and the momentum advection and pressure terms. 
    239262 
    240263A direct consequence of the pseudo-fourth order nature of the scheme is that it is not non-diffusive, 
    241264\ie\ the global variance of a tracer is not preserved using CEN4. 
    242 Furthermore, it must be used in conjunction with an explicit diffusion operator to produce a sensible solution. 
    243 As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction with an Asselin time-filter, 
    244 so $T$ in (\autoref{eq:TRA_adv_cen4}) is the \textit{now} tracer. 
     265Furthermore, it must be used in conjunction with an explicit diffusion operator to 
     266produce a sensible solution. 
     267As in CEN2 case, the time-stepping is performed using a leapfrog scheme in conjunction with 
     268an Asselin time-filter, so $T$ in (\autoref{eq:TRA_adv_cen4}) is the \textit{now} tracer. 
    245269 
    246270At a $T$-grid cell adjacent to a boundary (coastline, bottom and surface), 
     
    248272This hypothesis usually reduces the order of the scheme. 
    249273Here we choose to set the gradient of $T$ across the boundary to zero. 
    250 Alternative conditions can be specified, such as a reduction to a second order scheme for 
    251 these near boundary grid points. 
     274Alternative conditions can be specified, 
     275such as a reduction to a second order scheme for these near boundary grid points. 
    252276 
    253277%% ================================================================================================= 
     
    255279\label{subsec:TRA_adv_tvd} 
    256280 
    257 The Flux Corrected Transport schemes (FCT) is used when \np[=.true.]{ln_traadv_fct}{ln\_traadv\_fct}. 
    258 Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on horizontal (iso-level) and vertical direction by 
     281The \textbf{F}lux \textbf{C}orrected \textbf{T}ransport schemes (FCT) is used when 
     282\np[=.true.]{ln_traadv_fct}{ln\_traadv\_fct}. 
     283Its order ($2^{nd}$ or $4^{th}$) can be chosen independently on 
     284horizontal (iso-level) and vertical direction by 
    259285setting \np{nn_fct_h}{nn\_fct\_h} and \np{nn_fct_v}{nn\_fct\_v} to $2$ or $4$. 
    260286FCT implementation can be found in the \mdl{traadv\_fct} module. 
    261287 
    262 In FCT formulation, the tracer at velocity points is evaluated using a combination of an upstream and 
    263 a centred scheme. 
     288In FCT formulation, the tracer at velocity points is evaluated using 
     289a combination of an upstream and a centred scheme. 
    264290For example, in the $i$-direction : 
    265291\begin{equation} 
     
    270296                     T_{i + 1} & \text{if~} u_{i + 1/2} <    0 \\ 
    271297                     T_i       & \text{if~} u_{i + 1/2} \geq 0 \\ 
    272     \end{cases} 
    273     \\ 
     298    \end{cases} \\ 
    274299    \tau_u^{fct} &= \tau_u^{ups} + c_u \, \big( \tau_u^{cen} - \tau_u^{ups} \big) 
    275300  \end{split} 
     
    288313$\tau_u^{cen}$ is evaluated in (\autoref{eq:TRA_adv_fct}) using the \textit{now} tracer while 
    289314$\tau_u^{ups}$ is evaluated using the \textit{before} tracer. 
    290 In other words, the advective part of the scheme is time stepped with a leap-frog scheme 
    291 while a forward scheme is used for the diffusive part. 
     315In other words, the advective part of the scheme is time stepped with a leap-frog scheme while 
     316a forward scheme is used for the diffusive part. 
    292317 
    293318%% ================================================================================================= 
     
    295320\label{subsec:TRA_adv_mus} 
    296321 
    297 The Monotone Upstream Scheme for Conservative Laws (MUSCL) is used when \np[=.true.]{ln_traadv_mus}{ln\_traadv\_mus}. 
     322The \textbf{M}onotone \textbf{U}pstream \textbf{S}cheme for \textbf{C}onservative \textbf{L}aws 
     323(MUSCL) is used when \np[=.true.]{ln_traadv_mus}{ln\_traadv\_mus}. 
    298324MUSCL implementation can be found in the \mdl{traadv\_mus} module. 
    299325 
    300326MUSCL has been first implemented in \NEMO\ by \citet{levy.estublier.ea_GRL01}. 
    301 In its formulation, the tracer at velocity points is evaluated assuming a linear tracer variation between 
    302 two $T$-points (\autoref{fig:TRA_adv_scheme}). 
     327In its formulation, the tracer at velocity points is evaluated assuming 
     328a linear tracer variation between two $T$-points (\autoref{fig:TRA_adv_scheme}). 
    303329For example, in the $i$-direction : 
    304 \begin{equation} 
     330\[ 
    305331  % \label{eq:TRA_adv_mus} 
    306332  \tau_u^{mus} = \lt\{ 
    307333  \begin{split} 
    308                        \tau_i         &+ \frac{1}{2} \lt( 1 - \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
    309                        \widetilde{\partial_i         \tau} & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
    310                        \tau_{i + 1/2} &+ \frac{1}{2} \lt( 1 + \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
    311                        \widetilde{\partial_{i + 1/2} \tau} & \text{if~} u_{i + 1/2} <         0 
     334    \tau_i        &+ \frac{1}{2} \lt( 1 - \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
     335    \widetilde{\partial_i        \tau} & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
     336    \tau_{i + 1/2} &+ \frac{1}{2} \lt( 1 + \frac{u_{i + 1/2} \, \rdt}{e_{1u}} \rt) 
     337    \widetilde{\partial_{i + 1/2} \tau} & \text{if~} u_{i + 1/2} <         0 
    312338  \end{split} 
    313339                                                                                                      \rt. 
    314 \end{equation} 
    315 where $\widetilde{\partial_i \tau}$ is the slope of the tracer on which a limitation is imposed to 
    316 ensure the \textit{positive} character of the scheme. 
    317  
    318 The time stepping is performed using a forward scheme, that is the \textit{before} tracer field is used to 
    319 evaluate $\tau_u^{mus}$. 
     340\] 
     341where $\widetilde{\partial_i \tau}$ is the slope of the tracer on which 
     342a limitation is imposed to ensure the \textit{positive} character of the scheme. 
     343 
     344The time stepping is performed using a forward scheme, 
     345that is the \textit{before} tracer field is used to evaluate $\tau_u^{mus}$. 
    320346 
    321347For an ocean grid point adjacent to land and where the ocean velocity is directed toward land, 
    322348an upstream flux is used. 
    323349This choice ensure the \textit{positive} character of the scheme. 
    324 In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using upstream fluxes 
    325 (\np[=.true.]{ln_mus_ups}{ln\_mus\_ups}). 
     350In addition, fluxes round a grid-point where a runoff is applied can optionally be computed using 
     351upstream fluxes (\np[=.true.]{ln_mus_ups}{ln\_mus\_ups}). 
    326352 
    327353%% ================================================================================================= 
     
    329355\label{subsec:TRA_adv_ubs} 
    330356 
    331 The Upstream-Biased Scheme (UBS) is used when \np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs}. 
     357The \textbf{U}pstream-\textbf{B}iased \textbf{S}cheme (UBS) is used when 
     358\np[=.true.]{ln_traadv_ubs}{ln\_traadv\_ubs}. 
    332359UBS implementation can be found in the \mdl{traadv\_mus} module. 
    333360 
    334361The UBS scheme, often called UP3, is also known as the Cell Averaged QUICK scheme 
    335 (Quadratic Upstream Interpolation for Convective Kinematics). 
     362(\textbf{Q}uadratic \textbf{U}pstream \textbf{I}nterpolation for 
     363\textbf{C}onvective \textbf{K}inematics). 
    336364It is an upstream-biased third order scheme based on an upstream-biased parabolic interpolation. 
    337365For example, in the $i$-direction: 
     
    340368  \tau_u^{ubs} = \overline T ^{i + 1/2} - \frac{1}{6} 
    341369    \begin{cases} 
    342                                                       \tau"_i       & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
    343                                                       \tau"_{i + 1} & \text{if~} u_{i + 1/2} <         0 
     370      \tau"_i       & \text{if~} u_{i + 1/2} \geqslant 0 \\ 
     371      \tau"_{i + 1} & \text{if~} u_{i + 1/2} <         0 
    344372    \end{cases} 
    345   \quad 
    346   \text{where~} \tau"_i = \delta_i \lt[ \delta_{i + 1/2} [\tau] \rt] 
     373  \quad \text{where~} \tau"_i = \delta_i \lt[ \delta_{i + 1/2} [\tau] \rt] 
    347374\end{equation} 
    348375 
    349376This results in a dissipatively dominant (i.e. hyper-diffusive) truncation error 
    350377\citep{shchepetkin.mcwilliams_OM05}. 
    351 The overall performance of the advection scheme is similar to that reported in \cite{farrow.stevens_JPO95}. 
     378The overall performance of the advection scheme is similar to that reported in 
     379\cite{farrow.stevens_JPO95}. 
    352380It is a relatively good compromise between accuracy and smoothness. 
    353381Nevertheless the scheme is not \textit{positive}, meaning that false extrema are permitted, 
    354382but the amplitude of such are significantly reduced over the centred second or fourth order method. 
    355 Therefore it is not recommended that it should be applied to a passive tracer that requires positivity. 
     383Therefore it is not recommended that it should be applied to 
     384a passive tracer that requires positivity. 
    356385 
    357386The intrinsic diffusion of UBS makes its use risky in the vertical direction where 
    358387the control of artificial diapycnal fluxes is of paramount importance 
    359388\citep{shchepetkin.mcwilliams_OM05, demange_phd14}. 
    360 Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or a $4^th$ order COMPACT scheme 
    361 (\np[=2 or 4]{nn_ubs_v}{nn\_ubs\_v}). 
    362  
    363 For stability reasons (see \autoref{chap:TD}), the first term  in \autoref{eq:TRA_adv_ubs} 
    364 (which corresponds to a second order centred scheme) 
    365 is evaluated using the \textit{now} tracer (centred in time) while the second term 
    366 (which is the diffusive part of the scheme), 
     389Therefore the vertical flux is evaluated using either a $2^nd$ order FCT scheme or 
     390a $4^th$ order COMPACT scheme (\np[=2 or 4]{nn_ubs_v}{nn\_ubs\_v}). 
     391 
     392For stability reasons (see \autoref{chap:TD}), 
     393the first term  in \autoref{eq:TRA_adv_ubs} (which corresponds to a second order centred scheme) 
     394is evaluated using the \textit{now}    tracer (centred in time) while 
     395the second term (which is the diffusive part of the scheme), 
    367396is evaluated using the \textit{before} tracer (forward in time). 
    368 This choice is discussed by \citet{webb.de-cuevas.ea_JAOT98} in the context of the QUICK advection scheme. 
     397This choice is discussed by \citet{webb.de-cuevas.ea_JAOT98} in 
     398the context of the QUICK advection scheme. 
    369399UBS and QUICK schemes only differ by one coefficient. 
    370 Replacing 1/6 with 1/8 in \autoref{eq:TRA_adv_ubs} leads to the QUICK advection scheme \citep{webb.de-cuevas.ea_JAOT98}. 
     400Replacing 1/6 with 1/8 in \autoref{eq:TRA_adv_ubs} leads to the QUICK advection scheme 
     401\citep{webb.de-cuevas.ea_JAOT98}. 
    371402This option is not available through a namelist parameter, since the 1/6 coefficient is hard coded. 
    372 Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and obtain a QUICK scheme. 
     403Nevertheless it is quite easy to make the substitution in the \mdl{traadv\_ubs} module and 
     404obtain a QUICK scheme. 
    373405 
    374406Note that it is straightforward to rewrite \autoref{eq:TRA_adv_ubs} as follows: 
     
    389421Firstly, it clearly reveals that the UBS scheme is based on the fourth order scheme to which 
    390422an upstream-biased diffusion term is added. 
    391 Secondly, this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to 
    392 be evaluated at the \textit{now} time step using \autoref{eq:TRA_adv_ubs}. 
    393 Thirdly, the diffusion term is in fact a biharmonic operator with an eddy coefficient which 
    394 is simply proportional to the velocity: $A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$. 
    395 Note the current version of \NEMO\ uses the computationally more efficient formulation \autoref{eq:TRA_adv_ubs}. 
     423Secondly, 
     424this emphasises that the $4^{th}$ order part (as well as the $2^{nd}$ order part as stated above) has to be evaluated at the \textit{now} time step using \autoref{eq:TRA_adv_ubs}. 
     425Thirdly, the diffusion term is in fact a biharmonic operator with 
     426an eddy coefficient which is simply proportional to the velocity: 
     427$A_u^{lm} = \frac{1}{12} \, {e_{1u}}^3 \, |u|$. 
     428Note the current version of \NEMO\ uses the computationally more efficient formulation 
     429\autoref{eq:TRA_adv_ubs}. 
    396430 
    397431%% ================================================================================================= 
     
    399433\label{subsec:TRA_adv_qck} 
    400434 
    401 The Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms (QUICKEST) scheme 
    402 proposed by \citet{leonard_CMAME79} is used when \np[=.true.]{ln_traadv_qck}{ln\_traadv\_qck}. 
     435The \textbf{Q}uadratic \textbf{U}pstream \textbf{I}nterpolation for 
     436\textbf{C}onvective \textbf{K}inematics with \textbf{E}stimated \textbf{S}treaming \textbf{T}erms 
     437(QUICKEST) scheme proposed by \citet{leonard_CMAME79} is used when 
     438\np[=.true.]{ln_traadv_qck}{ln\_traadv\_qck}. 
    403439QUICKEST implementation can be found in the \mdl{traadv\_qck} module. 
    404440 
    405441QUICKEST is the third order Godunov scheme which is associated with the ULTIMATE QUICKEST limiter 
    406442\citep{leonard_CMAME91}. 
    407 It has been implemented in \NEMO\ by G. Reffray (MERCATOR-ocean) and can be found in the \mdl{traadv\_qck} module. 
     443It has been implemented in \NEMO\ by G. Reffray (Mercator Ocean) and 
     444can be found in the \mdl{traadv\_qck} module. 
    408445The resulting scheme is quite expensive but \textit{positive}. 
    409446It can be used on both active and passive tracers. 
     
    412449Therefore the vertical flux is evaluated using the CEN2 scheme. 
    413450This no longer guarantees the positivity of the scheme. 
    414 The use of FCT in the vertical direction (as for the UBS case) should be implemented to restore this property. 
    415  
    416 %%%gmcomment   :  Cross term are missing in the current implementation.... 
     451The use of FCT in the vertical direction (as for the UBS case) should be implemented to 
     452restore this property. 
     453 
     454\cmtgm{Cross term are missing in the current implementation....} 
    417455 
    418456%% ================================================================================================= 
     
    428466Options are defined through the \nam{tra_ldf}{tra\_ldf} namelist variables. 
    429467They are regrouped in four items, allowing to specify 
    430 $(i)$   the type of operator used (none, laplacian, bilaplacian), 
    431 $(ii)$  the direction along which the operator acts (iso-level, horizontal, iso-neutral), 
    432 $(iii)$ some specific options related to the rotated operators (\ie\ non-iso-level operator), and 
    433 $(iv)$  the specification of eddy diffusivity coefficient (either constant or variable in space and time). 
    434 Item $(iv)$ will be described in \autoref{chap:LDF}. 
     468\begin{enumerate*}[label=(\textit{\roman*})] 
     469\item the type of operator used (none, laplacian, bilaplacian), 
     470\item the direction along which the operator acts (iso-level, horizontal, iso-neutral), 
     471\item some specific options related to the rotated operators (\ie\ non-iso-level operator), and 
     472\item the specification of eddy diffusivity coefficient 
     473  (either constant or variable in space and time). 
     474\end{enumerate*} 
     475Item (iv) will be described in \autoref{chap:LDF}. 
    435476The direction along which the operators act is defined through the slope between 
    436477this direction and the iso-level surfaces. 
     
    440481\ie\ the tracers appearing in its expression are the \textit{before} tracers in time, 
    441482except for the pure vertical component that appears when a rotation tensor is used. 
    442 This latter component is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}). 
    443 When \np[=.true.]{ln_traldf_msc}{ln\_traldf\_msc}, a Method of Stabilizing Correction is used in which 
    444 the pure vertical component is split into an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}. 
     483This latter component is solved implicitly together with the vertical diffusion term 
     484(see \autoref{chap:TD}). 
     485When \np[=.true.]{ln_traldf_msc}{ln\_traldf\_msc}, 
     486a Method of Stabilizing Correction is used in which the pure vertical component is split into 
     487an explicit and an implicit part \citep{lemarie.debreu.ea_OM12}. 
    445488 
    446489%% ================================================================================================= 
     
    451494 
    452495\begin{description} 
    453 \item [{\np[=.true.]{ln_traldf_OFF}{ln\_traldf\_OFF}}] no operator selected, the lateral diffusive tendency will not be applied to the tracer equation. 
    454   This option can be used when the selected advection scheme is diffusive enough (MUSCL scheme for example). 
     496\item [{\np[=.true.]{ln_traldf_OFF}{ln\_traldf\_OFF}}] no operator selected, 
     497  the lateral diffusive tendency will not be applied to the tracer equation. 
     498  This option can be used when the selected advection scheme is diffusive enough 
     499  (MUSCL scheme for example). 
    455500\item [{\np[=.true.]{ln_traldf_lap}{ln\_traldf\_lap}}] a laplacian operator is selected. 
    456   This harmonic operator takes the following expression:  $\mathcal{L}(T) = \nabla \cdot A_{ht} \; \nabla T $, 
     501  This harmonic operator takes the following expression: 
     502  $\mathcal{L}(T) = \nabla \cdot A_{ht} \; \nabla T $, 
    457503  where the gradient operates along the selected direction (see \autoref{subsec:TRA_ldf_dir}), 
    458504  and $A_{ht}$ is the eddy diffusivity coefficient expressed in $m^2/s$ (see \autoref{chap:LDF}). 
     
    461507  $\mathcal{B} = - \mathcal{L}(\mathcal{L}(T)) = - \nabla \cdot b \nabla (\nabla \cdot b \nabla T)$ 
    462508  where the gradient operats along the selected direction, 
    463   and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ (see \autoref{chap:LDF}). 
     509  and $b^2 = B_{ht}$ is the eddy diffusivity coefficient expressed in $m^4/s$ 
     510  (see \autoref{chap:LDF}). 
    464511  In the code, the bilaplacian operator is obtained by calling the laplacian twice. 
    465512\end{description} 
     
    469516minimizing the impact on the larger scale features. 
    470517The main difference between the two operators is the scale selectiveness. 
    471 The bilaplacian damping time (\ie\ its spin down time) scales like $\lambda^{-4}$ for 
    472 disturbances of wavelength $\lambda$ (so that short waves damped more rapidelly than long ones), 
     518The bilaplacian damping time (\ie\ its spin down time) scales like 
     519$\lambda^{-4}$ for disturbances of wavelength $\lambda$ 
     520(so that short waves damped more rapidelly than long ones), 
    473521whereas the laplacian damping time scales only like $\lambda^{-2}$. 
    474522 
     
    479527The choice of a direction of action determines the form of operator used. 
    480528The operator is a simple (re-entrant) laplacian acting in the (\textbf{i},\textbf{j}) plane when 
    481 iso-level option is used (\np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev}) or 
    482 when a horizontal (\ie\ geopotential) operator is demanded in \textit{z}-coordinate 
     529iso-level option is used (\np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev}) or when 
     530a horizontal (\ie\ geopotential) operator is demanded in \textit{z}-coordinate 
    483531(\np{ln_traldf_hor}{ln\_traldf\_hor} and \np[=.true.]{ln_zco}{ln\_zco}). 
    484532The associated code can be found in the \mdl{traldf\_lap\_blp} module. 
     
    489537see \mdl{traldf\_iso} or \mdl{traldf\_triad} module, resp.), or 
    490538when a horizontal (\ie\ geopotential) operator is demanded in \textit{s}-coordinate 
    491 (\np{ln_traldf_hor}{ln\_traldf\_hor} and \np{ln_sco}{ln\_sco} = \forcode{.true.}) 
    492 \footnote{In this case, the standard iso-neutral operator will be automatically selected}. 
     539(\np{ln_traldf_hor}{ln\_traldf\_hor} and \np{ln_sco}{ln\_sco} = \forcode{.true.}) \footnote{ 
     540  In this case, the standard iso-neutral operator will be automatically selected}. 
    493541In that case, a rotation is applied to the gradient(s) that appears in the operator so that 
    494542diffusive fluxes acts on the three spatial direction. 
     
    511559first (and third in bilaplacian case) horizontal tracer derivative are masked. 
    512560It is implemented in the \rou{tra\_ldf\_lap} subroutine found in the \mdl{traldf\_lap\_blp} module. 
    513 The module also contains \rou{tra\_ldf\_blp}, the subroutine calling twice \rou{tra\_ldf\_lap} in order to 
     561The module also contains \rou{tra\_ldf\_blp}, 
     562the subroutine calling twice \rou{tra\_ldf\_lap} in order to 
    514563compute the iso-level bilaplacian operator. 
    515564 
    516565It is a \textit{horizontal} operator (\ie acting along geopotential surfaces) in 
    517 the $z$-coordinate with or without partial steps, but is simply an iso-level operator in the $s$-coordinate. 
    518 It is thus used when, in addition to \np{ln_traldf_lap}{ln\_traldf\_lap} or \np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp}, 
    519 we have \np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev} or \np{ln_traldf_hor}{ln\_traldf\_hor}~=~\np[=.true.]{ln_zco}{ln\_zco}. 
     566the $z$-coordinate with or without partial steps, 
     567but is simply an iso-level operator in the $s$-coordinate. 
     568It is thus used when, 
     569in addition to \np{ln_traldf_lap}{ln\_traldf\_lap} or \np[=.true.]{ln_traldf_blp}{ln\_traldf\_blp}, 
     570we have \np[=.true.]{ln_traldf_lev}{ln\_traldf\_lev} or 
     571\np[=]{ln_traldf_hor}{ln\_traldf\_hor}\np[=.true.]{ln_zco}{ln\_zco}. 
    520572In both cases, it significantly contributes to diapycnal mixing. 
    521573It is therefore never recommended, even when using it in the bilaplacian case. 
     
    523575Note that in the partial step $z$-coordinate (\np[=.true.]{ln_zps}{ln\_zps}), 
    524576tracers in horizontally adjacent cells are located at different depths in the vicinity of the bottom. 
    525 In this case, horizontal derivatives in (\autoref{eq:TRA_ldf_lap}) at the bottom level require a specific treatment. 
     577In this case, 
     578horizontal derivatives in (\autoref{eq:TRA_ldf_lap}) at the bottom level require a specific treatment. 
    526579They are calculated in the \mdl{zpshde} module, described in \autoref{sec:TRA_zpshde}. 
    527580 
     
    533586\subsubsection[Standard rotated (bi-)laplacian operator (\textit{traldf\_iso.F90})]{Standard rotated (bi-)laplacian operator (\protect\mdl{traldf\_iso})} 
    534587\label{subsec:TRA_ldf_iso} 
     588 
    535589The general form of the second order lateral tracer subgrid scale physics (\autoref{eq:MB_zdf}) 
    536 takes the following semi -discrete space form in $z$- and $s$-coordinates: 
     590takes the following semi-discrete space form in $z$- and $s$-coordinates: 
    537591\begin{equation} 
    538592  \label{eq:TRA_ldf_iso} 
     
    554608or both \np[=.true.]{ln_traldf_hor}{ln\_traldf\_hor} and \np[=.true.]{ln_zco}{ln\_zco}. 
    555609The way these slopes are evaluated is given in \autoref{sec:LDF_slp}. 
    556 At the surface, bottom and lateral boundaries, the turbulent fluxes of heat and salt are set to zero using 
    557 the mask technique (see \autoref{sec:LBC_coast}). 
     610At the surface, bottom and lateral boundaries, 
     611the turbulent fluxes of heat and salt are set to zero using the mask technique 
     612(see \autoref{sec:LBC_coast}). 
    558613 
    559614The operator in \autoref{eq:TRA_ldf_iso} involves both lateral and vertical derivatives. 
    560 For numerical stability, the vertical second derivative must be solved using the same implicit time scheme as that 
    561 used in the vertical physics (see \autoref{sec:TRA_zdf}). 
     615For numerical stability, the vertical second derivative must be solved using 
     616the same implicit time scheme as that used in the vertical physics (see \autoref{sec:TRA_zdf}). 
    562617For computer efficiency reasons, this term is not computed in the \mdl{traldf\_iso} module, 
    563618but in the \mdl{trazdf} module where, if iso-neutral mixing is used, 
    564 the vertical mixing coefficient is simply increased by $\frac{e_{1w} e_{2w}}{e_{3w}}(r_{1w}^2 + r_{2w}^2)$. 
     619the vertical mixing coefficient is simply increased by 
     620$\frac{e_{1w} e_{2w}}{e_{3w}}(r_{1w}^2 + r_{2w}^2)$. 
    565621 
    566622This formulation conserves the tracer but does not ensure the decrease of the tracer variance. 
    567 Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to run safely without 
    568 any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}. 
     623Nevertheless the treatment performed on the slopes (see \autoref{chap:LDF}) allows the model to 
     624run safely without any additional background horizontal diffusion \citep{guilyardi.madec.ea_CD01}. 
    569625 
    570626Note that in the partial step $z$-coordinate (\np[=.true.]{ln_zps}{ln\_zps}), 
    571 the horizontal derivatives at the bottom level in \autoref{eq:TRA_ldf_iso} require a specific treatment. 
     627the horizontal derivatives at the bottom level in \autoref{eq:TRA_ldf_iso} require 
     628a specific treatment. 
    572629They are calculated in module zpshde, described in \autoref{sec:TRA_zpshde}. 
    573630 
     
    576633\label{subsec:TRA_ldf_triad} 
    577634 
    578 An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which ensures tracer variance decreases 
    579 is also available in \NEMO\ (\np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}). 
     635An alternative scheme developed by \cite{griffies.gnanadesikan.ea_JPO98} which 
     636ensures tracer variance decreases is also available in \NEMO\ 
     637(\np[=.true.]{ln_traldf_triad}{ln\_traldf\_triad}). 
    580638A complete description of the algorithm is given in \autoref{apdx:TRIADS}. 
    581639 
    582 The lateral fourth order bilaplacian operator on tracers is obtained by applying (\autoref{eq:TRA_ldf_lap}) twice. 
     640The lateral fourth order bilaplacian operator on tracers is obtained by 
     641applying (\autoref{eq:TRA_ldf_lap}) twice. 
    583642The operator requires an additional assumption on boundary conditions: 
    584643both first and third derivative terms normal to the coast are set to zero. 
    585644 
    586 The lateral fourth order operator formulation on tracers is obtained by applying (\autoref{eq:TRA_ldf_iso}) twice. 
     645The lateral fourth order operator formulation on tracers is obtained by 
     646applying (\autoref{eq:TRA_ldf_iso}) twice. 
    587647It requires an additional assumption on boundary conditions: 
    588648first and third derivative terms normal to the coast, 
     
    593653\label{subsec:TRA_ldf_options} 
    594654 
    595 \begin{itemize} 
    596 \item \np{ln_traldf_msc}{ln\_traldf\_msc} = Method of Stabilizing Correction (both operators) 
    597 \item \np{rn_slpmax}{rn\_slpmax} = slope limit (both operators) 
    598 \item \np{ln_triad_iso}{ln\_triad\_iso} = pure horizontal mixing in ML (triad only) 
    599 \item \np{rn_sw_triad}{rn\_sw\_triad} $= 1$ switching triad; $= 0$ all 4 triads used (triad only) 
    600 \item \np{ln_botmix_triad}{ln\_botmix\_triad} = lateral mixing on bottom (triad only) 
    601 \end{itemize} 
     655\begin{labeling}{{\np{ln_botmix_triad}{ln\_botmix\_triad}}} 
     656\item [{\np{ln_traldf_msc}{ln\_traldf\_msc}    }] Method of Stabilizing Correction (both operators) 
     657\item [{\np{rn_slpmax}{rn\_slpmax}             }] Slope limit (both operators) 
     658\item [{\np{ln_triad_iso}{ln\_triad\_iso}      }] Pure horizontal mixing in ML (triad only) 
     659\item [{\np{rn_sw_triad}{rn\_sw\_triad}        }] \forcode{=1} switching triad; 
     660  \forcode{= 0} all 4 triads used (triad only) 
     661\item [{\np{ln_botmix_triad}{ln\_botmix\_triad}}] Lateral mixing on bottom (triad only) 
     662\end{labeling} 
    602663 
    603664%% ================================================================================================= 
     
    606667 
    607668Options are defined through the \nam{zdf}{zdf} namelist variables. 
    608 The formulation of the vertical subgrid scale tracer physics is the same for all the vertical coordinates, 
    609 and is based on a laplacian operator. 
    610 The vertical diffusion operator given by (\autoref{eq:MB_zdf}) takes the following semi -discrete space form: 
    611 \begin{gather*} 
     669The formulation of the vertical subgrid scale tracer physics is the same for 
     670all the vertical coordinates, and is based on a laplacian operator. 
     671The vertical diffusion operator given by (\autoref{eq:MB_zdf}) takes 
     672the following semi-discrete space form: 
     673\[ 
    612674  % \label{eq:TRA_zdf} 
    613     D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \\ 
    614     D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt] 
    615 \end{gather*} 
    616 where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on temperature and salinity, 
    617 respectively. 
     675  D^{vT}_T = \frac{1}{e_{3t}} \, \delta_k \lt[ \, \frac{A^{vT}_w}{e_{3w}} \delta_{k + 1/2}[T] \, \rt] \quad 
     676  D^{vS}_T = \frac{1}{e_{3t}} \; \delta_k \lt[ \, \frac{A^{vS}_w}{e_{3w}} \delta_{k + 1/2}[S] \, \rt] 
     677\] 
     678where $A_w^{vT}$ and $A_w^{vS}$ are the vertical eddy diffusivity coefficients on 
     679temperature and salinity, respectively. 
    618680Generally, $A_w^{vT} = A_w^{vS}$ except when double diffusive mixing is parameterised 
    619681(\ie\ \np[=.true.]{ln_zdfddm}{ln\_zdfddm},). 
    620682The way these coefficients are evaluated is given in \autoref{chap:ZDF} (ZDF). 
    621 Furthermore, when iso-neutral mixing is used, both mixing coefficients are increased by 
    622 $\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to account for the vertical second derivative of 
    623 \autoref{eq:TRA_ldf_iso}. 
     683Furthermore, when iso-neutral mixing is used, 
     684both mixing coefficients are increased by $\frac{e_{1w} e_{2w}}{e_{3w} }({r_{1w}^2 + r_{2w}^2})$ to 
     685account for the vertical second derivative of \autoref{eq:TRA_ldf_iso}. 
    624686 
    625687At the surface and bottom boundaries, the turbulent fluxes of heat and salt must be specified. 
     
    628690a geothermal flux forcing is prescribed as a bottom boundary condition (see \autoref{subsec:TRA_bbc}). 
    629691 
    630 The large eddy coefficient found in the mixed layer together with high vertical resolution implies that 
    631 there would be too restrictive constraint on the time step if we use explicit time stepping. 
     692The large eddy coefficient found in the mixed layer together with high vertical resolution implies 
     693that there would be too restrictive constraint on the time step if we use explicit time stepping. 
    632694Therefore an implicit time stepping is preferred for the vertical diffusion since 
    633695it overcomes the stability constraint. 
     
    648710 
    649711Due to interactions and mass exchange of water ($F_{mass}$) with other Earth system components 
    650 (\ie\ atmosphere, sea-ice, land), the change in the heat and salt content of the surface layer of the ocean is due 
    651 both to the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) and 
     712(\ie\ atmosphere, sea-ice, land), 
     713the change in the heat and salt content of the surface layer of the ocean is due both to 
     714the heat and salt fluxes crossing the sea surface (not linked with $F_{mass}$) and 
    652715to the heat and salt content of the mass exchange. 
    653716They are both included directly in $Q_{ns}$, the surface heat flux, 
    654717and $F_{salt}$, the surface salt flux (see \autoref{chap:SBC} for further details). 
    655 By doing this, the forcing formulation is the same for any tracer (including temperature and salinity). 
    656  
    657 The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields (used on tracers): 
    658  
    659 \begin{itemize} 
    660 \item $Q_{ns}$, the non-solar part of the net surface heat flux that crosses the sea surface 
    661   (\ie\ the difference between the total surface heat flux and the fraction of the short wave flux that 
    662   penetrates into the water column, see \autoref{subsec:TRA_qsr}) 
     718By doing this, the forcing formulation is the same for any tracer 
     719(including temperature and salinity). 
     720 
     721The surface module (\mdl{sbcmod}, see \autoref{chap:SBC}) provides the following forcing fields 
     722(used on tracers): 
     723 
     724\begin{labeling}{\textit{fwfisf}} 
     725\item [$Q_{ns}$] The non-solar part of the net surface heat flux that crosses the sea surface 
     726  (\ie\ the difference between the total surface heat flux and 
     727  the fraction of the short wave flux that penetrates into the water column, 
     728  see \autoref{subsec:TRA_qsr}) 
    663729  plus the heat content associated with of the mass exchange with the atmosphere and lands. 
    664 \item $\textit{sfx}$, the salt flux resulting from ice-ocean mass exchange (freezing, melting, ridging...) 
    665 \item \textit{emp}, the mass flux exchanged with the atmosphere (evaporation minus precipitation) and 
     730\item [\textit{sfx}] The salt flux resulting from ice-ocean mass exchange 
     731  (freezing, melting, ridging...) 
     732\item [\textit{emp}] The mass flux exchanged with the atmosphere (evaporation minus precipitation) and 
    666733  possibly with the sea-ice and ice-shelves. 
    667 \item \textit{rnf}, the mass flux associated with runoff 
     734\item [\textit{rnf}] The mass flux associated with runoff 
    668735  (see \autoref{sec:SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
    669 \item \textit{fwfisf}, the mass flux associated with ice shelf melt, 
     736\item [\textit{fwfisf}] The mass flux associated with ice shelf melt, 
    670737  (see \autoref{sec:SBC_isf} for further details on how the ice shelf melt is computed and applied). 
    671 \end{itemize} 
     738\end{labeling} 
    672739 
    673740The surface boundary condition on temperature and salinity is applied as follows: 
    674741\begin{equation} 
    675742  \label{eq:TRA_sbc} 
    676   \begin{alignedat}{2} 
    677     F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{Q_{ns}      }^t \\ 
    678     F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} &\overline{\textit{sfx}}^t 
    679   \end{alignedat} 
     743    F^T = \frac{1}{C_p} \frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} \overline{Q_{ns}      }^t \qquad 
     744    F^S =               \frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} \overline{\textit{sfx}}^t 
    680745\end{equation} 
    681746where $\overline x^t$ means that $x$ is averaged over two consecutive time steps 
     
    683748Such time averaging prevents the divergence of odd and even time step (see \autoref{chap:TD}). 
    684749 
    685 In the linear free surface case (\np[=.true.]{ln_linssh}{ln\_linssh}), an additional term has to be added on 
    686 both temperature and salinity. 
    687 On temperature, this term remove the heat content associated with mass exchange that has been added to $Q_{ns}$. 
    688 On salinity, this term mimics the concentration/dilution effect that would have resulted from a change in 
    689 the volume of the first level. 
     750In the linear free surface case (\np[=.true.]{ln_linssh}{ln\_linssh}), 
     751an additional term has to be added on both temperature and salinity. 
     752On temperature, this term remove the heat content associated with 
     753mass exchange that has been added to $Q_{ns}$. 
     754On salinity, this term mimics the concentration/dilution effect that would have resulted from 
     755a change in the volume of the first level. 
    690756The resulting surface boundary condition is applied as follows: 
    691757\begin{equation} 
    692758  \label{eq:TRA_sbc_lin} 
    693   \begin{alignedat}{2} 
    694     F^T &= \frac{1}{C_p} &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
    695           &\overline{(Q_{ns}       - C_p \, \textit{emp} \lt. T \rt|_{k = 1})}^t \\ 
    696     F^S &=               &\frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
    697           &\overline{(\textit{sfx} -        \textit{emp} \lt. S \rt|_{k = 1})}^t 
    698   \end{alignedat} 
    699 \end{equation} 
    700 Note that an exact conservation of heat and salt content is only achieved with non-linear free surface. 
     759    F^T = \frac{1}{C_p} \frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
     760          \overline{(Q_{ns}       - C_p \, \textit{emp} \lt. T \rt|_{k = 1})}^t \qquad 
     761    F^S =               \frac{1}{\rho_o \lt. e_{3t} \rt|_{k = 1}} 
     762          \overline{(\textit{sfx} -        \textit{emp} \lt. S \rt|_{k = 1})}^t 
     763\end{equation} 
     764Note that an exact conservation of heat and salt content is only achieved with 
     765non-linear free surface. 
    701766In the linear free surface case, there is a small imbalance. 
    702 The imbalance is larger than the imbalance associated with the Asselin time filter \citep{leclair.madec_OM09}. 
    703 This is the reason why the modified filter is not applied in the linear free surface case (see \autoref{chap:TD}). 
     767The imbalance is larger than the imbalance associated with the Asselin time filter 
     768\citep{leclair.madec_OM09}. 
     769This is the reason why the modified filter is not applied in the linear free surface case 
     770(see \autoref{chap:TD}). 
    704771 
    705772%% ================================================================================================= 
     
    716783When the penetrative solar radiation option is used (\np[=.true.]{ln_traqsr}{ln\_traqsr}), 
    717784the solar radiation penetrates the top few tens of meters of the ocean. 
    718 If it is not used (\np[=.false.]{ln_traqsr}{ln\_traqsr}) all the heat flux is absorbed in the first ocean level. 
    719 Thus, in the former case a term is added to the time evolution equation of temperature \autoref{eq:MB_PE_tra_T} and 
    720 the surface boundary condition is modified to take into account only the non-penetrative part of the surface 
    721 heat flux: 
     785If it is not used (\np[=.false.]{ln_traqsr}{ln\_traqsr}) all the heat flux is absorbed in 
     786the first ocean level. 
     787Thus, in the former case a term is added to the time evolution equation of temperature 
     788\autoref{eq:MB_PE_tra_T} and the surface boundary condition is modified to 
     789take into account only the non-penetrative part of the surface heat flux: 
    722790\begin{equation} 
    723791  \label{eq:TRA_PE_qsr} 
     
    736804 
    737805The shortwave radiation, $Q_{sr}$, consists of energy distributed across a wide spectral range. 
    738 The ocean is strongly absorbing for wavelengths longer than 700~nm and these wavelengths contribute to 
    739 heating the upper few tens of centimetres. 
    740 The fraction of $Q_{sr}$ that resides in these almost non-penetrative wavebands, $R$, is $\sim 58\%$ 
     806The ocean is strongly absorbing for wavelengths longer than 700 $nm$ and 
     807these wavelengths contribute to heat the upper few tens of centimetres. 
     808The fraction of $Q_{sr}$ that resides in these almost non-penetrative wavebands, $R$, is $\sim$ 58\% 
    741809(specified through namelist parameter \np{rn_abs}{rn\_abs}). 
    742 It is assumed to penetrate the ocean with a decreasing exponential profile, with an e-folding depth scale, $\xi_0$, 
    743 of a few tens of centimetres (typically $\xi_0 = 0.35~m$ set as \np{rn_si0}{rn\_si0} in the \nam{tra_qsr}{tra\_qsr} namelist). 
    744 For shorter wavelengths (400-700~nm), the ocean is more transparent, and solar energy propagates to 
    745 larger depths where it contributes to local heating. 
    746 The way this second part of the solar energy penetrates into the ocean depends on which formulation is chosen. 
     810It is assumed to penetrate the ocean with a decreasing exponential profile, 
     811with an e-folding depth scale, $\xi_0$, of a few tens of centimetres 
     812(typically $\xi_0 = 0.35~m$ set as \np{rn_si0}{rn\_si0} in the \nam{tra_qsr}{tra\_qsr} namelist). 
     813For shorter wavelengths (400-700 $nm$), the ocean is more transparent, 
     814and solar energy propagates to larger depths where it contributes to local heating. 
     815The way this second part of the solar energy penetrates into 
     816the ocean depends on which formulation is chosen. 
    747817In the simple 2-waveband light penetration scheme (\np[=.true.]{ln_qsr_2bd}{ln\_qsr\_2bd}) 
    748818a chlorophyll-independent monochromatic formulation is chosen for the shorter wavelengths, 
     
    754824where $\xi_1$ is the second extinction length scale associated with the shorter wavelengths. 
    755825It is usually chosen to be 23~m by setting the \np{rn_si0}{rn\_si0} namelist parameter. 
    756 The set of default values ($\xi_0, \xi_1, R$) corresponds to a Type I water in Jerlov's (1968) classification 
    757 (oligotrophic waters). 
     826The set of default values ($\xi_0, \xi_1, R$) corresponds to 
     827a Type I water in Jerlov's (1968) classification (oligotrophic waters). 
    758828 
    759829Such assumptions have been shown to provide a very crude and simplistic representation of 
     
    763833a 61 waveband formulation. 
    764834Unfortunately, such a model is very computationally expensive. 
    765 Thus, \cite{lengaigne.menkes.ea_CD07} have constructed a simplified version of this formulation in which 
    766 visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700nm). 
    767 For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to the coefficients computed from 
    768 the full spectral model of \cite{morel_JGR88} (as modified by \cite{morel.maritorena_JGR01}), 
    769 assuming the same power-law relationship. 
    770 As shown in \autoref{fig:TRA_qsr_irradiance}, this formulation, called RGB (Red-Green-Blue), 
     835Thus, \cite{lengaigne.menkes.ea_CD07} have constructed a simplified version of 
     836this formulation in which visible light is split into three wavebands: 
     837blue (400-500 $nm$), green (500-600 $nm$) and red (600-700 $nm$). 
     838For each wave-band, the chlorophyll-dependent attenuation coefficient is fitted to 
     839the coefficients computed from the full spectral model of \cite{morel_JGR88} 
     840(as modified by \cite{morel.maritorena_JGR01}), assuming the same power-law relationship. 
     841As shown in \autoref{fig:TRA_qsr_irradiance}, this formulation, 
     842called RGB (\textbf{R}ed-\textbf{G}reen-\textbf{B}lue), 
    771843reproduces quite closely the light penetration profiles predicted by the full spectal model, 
    772844but with much greater computational efficiency. 
     
    774846 
    775847The RGB formulation is used when \np[=.true.]{ln_qsr_rgb}{ln\_qsr\_rgb}. 
    776 The RGB attenuation coefficients (\ie\ the inverses of the extinction length scales) are tabulated over 
    777 61 nonuniform chlorophyll classes ranging from 0.01 to 10 g.Chl/L 
     848The RGB attenuation coefficients (\ie\ the inverses of the extinction length scales) are 
     849tabulated over 61 nonuniform chlorophyll classes ranging from 0.01 to 10 $g.Chl/L$ 
    778850(see the routine \rou{trc\_oce\_rgb} in \mdl{trc\_oce} module). 
    779851Four types of chlorophyll can be chosen in the RGB formulation: 
    780852 
    781853\begin{description} 
    782 \item [{\np[=0]{nn_chldta}{nn\_chldta}}] a constant 0.05 g.Chl/L value everywhere ; 
    783 \item [{\np[=1]{nn_chldta}{nn\_chldta}}] an observed time varying chlorophyll deduced from satellite surface ocean color measurement spread uniformly in the vertical direction; 
    784 \item [{\np[=2]{nn_chldta}{nn\_chldta}}] same as previous case except that a vertical profile of chlorophyl is used. 
    785   Following \cite{morel.berthon_LO89}, the profile is computed from the local surface chlorophyll value; 
    786 \item [{\np[=.true.]{ln_qsr_bio}{ln\_qsr\_bio}}] simulated time varying chlorophyll by TOP biogeochemical model. 
    787   In this case, the RGB formulation is used to calculate both the phytoplankton light limitation in 
    788   PISCES and the oceanic heating rate. 
     854\item [{\np[=0]{nn_chldta}{nn\_chldta}}] a constant 0.05 $g.Chl/L$ value everywhere; 
     855\item [{\np[=1]{nn_chldta}{nn\_chldta}}] an observed time varying chlorophyll deduced from 
     856  satellite surface ocean color measurement spread uniformly in the vertical direction; 
     857\item [{\np[=2]{nn_chldta}{nn\_chldta}}] same as previous case except that 
     858  a vertical profile of chlorophyl is used. 
     859  Following \cite{morel.berthon_LO89}, 
     860  the profile is computed from the local surface chlorophyll value; 
     861\item [{\np[=.true.]{ln_qsr_bio}{ln\_qsr\_bio}}] simulated time varying chlorophyll by 
     862  \TOP\ biogeochemical model. 
     863  In this case, the RGB formulation is used to calculate both 
     864  the phytoplankton light limitation in \PISCES\ and the oceanic heating rate. 
    789865\end{description} 
    790866 
     
    797873(\ie\ it is less than the computer precision) is computed once, 
    798874and the trend associated with the penetration of the solar radiation is only added down to that level. 
    799 Finally, note that when the ocean is shallow ($<$ 200~m), part of the solar radiation can reach the ocean floor. 
     875Finally, note that when the ocean is shallow ($<$ 200~m), 
     876part of the solar radiation can reach the ocean floor. 
    800877In this case, we have chosen that all remaining radiation is absorbed in the last ocean level 
    801878(\ie\ $I$ is masked). 
    802879 
    803 \begin{figure}[!t] 
     880\begin{figure} 
    804881  \centering 
    805   \includegraphics[width=0.66\textwidth]{Fig_TRA_Irradiance} 
     882  \includegraphics[width=0.66\textwidth]{TRA_Irradiance} 
    806883  \caption[Penetration profile of the downward solar irradiance calculated by four models]{ 
    807884    Penetration profile of the downward solar irradiance calculated by four models. 
     
    810887    4 waveband RGB formulation (red), 
    811888    61 waveband Morel (1988) formulation (black) for a chlorophyll concentration of 
    812     (a) Chl=0.05 mg/m$^3$ and (b) Chl=0.5 mg/m$^3$. 
     889    (a) Chl=0.05 $mg/m^3$ and (b) Chl=0.5 $mg/m^3$. 
    813890    From \citet{lengaigne.menkes.ea_CD07}.} 
    814891  \label{fig:TRA_qsr_irradiance} 
     
    824901  \label{lst:nambbc} 
    825902\end{listing} 
    826 \begin{figure}[!t] 
     903 
     904\begin{figure} 
    827905  \centering 
    828   \includegraphics[width=0.66\textwidth]{Fig_TRA_geoth} 
     906  \includegraphics[width=0.66\textwidth]{TRA_geoth} 
    829907  \caption[Geothermal heat flux]{ 
    830908    Geothermal Heat flux (in $mW.m^{-2}$) used by \cite{emile-geay.madec_OS09}. 
     
    836914\ie\ a no flux boundary condition is applied on active tracers at the bottom. 
    837915This is the default option in \NEMO, and it is implemented using the masking technique. 
    838 However, there is a non-zero heat flux across the seafloor that is associated with solid earth cooling. 
    839 This flux is weak compared to surface fluxes (a mean global value of $\sim 0.1 \, W/m^2$ \citep{stein.stein_N92}), 
     916However, there is a non-zero heat flux across the seafloor that 
     917is associated with solid earth cooling. 
     918This flux is weak compared to surface fluxes 
     919(a mean global value of $\sim 0.1 \, W/m^2$ \citep{stein.stein_N92}), 
    840920but it warms systematically the ocean and acts on the densest water masses. 
    841921Taking this flux into account in a global ocean model increases the deepest overturning cell 
    842 (\ie\ the one associated with the Antarctic Bottom Water) by a few Sverdrups \citep{emile-geay.madec_OS09}. 
     922(\ie\ the one associated with the Antarctic Bottom Water) by 
     923a few Sverdrups \citep{emile-geay.madec_OS09}. 
    843924 
    844925Options are defined through the \nam{bbc}{bbc} namelist variables. 
    845 The presence of geothermal heating is controlled by setting the namelist parameter \np{ln_trabbc}{ln\_trabbc} to true. 
    846 Then, when \np{nn_geoflx}{nn\_geoflx} is set to 1, a constant geothermal heating is introduced whose value is given by 
    847 the \np{rn_geoflx_cst}{rn\_geoflx\_cst}, which is also a namelist parameter. 
    848 When \np{nn_geoflx}{nn\_geoflx} is set to 2, a spatially varying geothermal heat flux is introduced which is provided in 
    849 the \ifile{geothermal\_heating} NetCDF file (\autoref{fig:TRA_geothermal}) \citep{emile-geay.madec_OS09}. 
     926The presence of geothermal heating is controlled by 
     927setting the namelist parameter \np{ln_trabbc}{ln\_trabbc} to true. 
     928Then, when \np{nn_geoflx}{nn\_geoflx} is set to 1, a constant geothermal heating is introduced whose 
     929value is given by the \np{rn_geoflx_cst}{rn\_geoflx\_cst}, which is also a namelist parameter. 
     930When \np{nn_geoflx}{nn\_geoflx} is set to 2, 
     931a spatially varying geothermal heat flux is introduced which is provided in 
     932the \ifile{geothermal\_heating} NetCDF file 
     933(\autoref{fig:TRA_geothermal}) \citep{emile-geay.madec_OS09}. 
    850934 
    851935%% ================================================================================================= 
     
    865949where dense water formed in marginal seas flows into a basin filled with less dense water, 
    866950or along the continental slope when dense water masses are formed on a continental shelf. 
    867 The amount of entrainment that occurs in these gravity plumes is critical in determining the density and 
    868 volume flux of the densest waters of the ocean, such as Antarctic Bottom Water, or North Atlantic Deep Water. 
     951The amount of entrainment that occurs in these gravity plumes is critical in 
     952determining the density and volume flux of the densest waters of the ocean, 
     953such as Antarctic Bottom Water, or North Atlantic Deep Water. 
    869954$z$-coordinate models tend to overestimate the entrainment, 
    870 because the gravity flow is mixed vertically by convection as it goes ''downstairs'' following the step topography, 
     955because the gravity flow is mixed vertically by convection as 
     956it goes ''downstairs'' following the step topography, 
    871957sometimes over a thickness much larger than the thickness of the observed gravity plume. 
    872 A similar problem occurs in the $s$-coordinate when the thickness of the bottom level varies rapidly downstream of 
    873 a sill \citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved. 
    874  
    875 The idea of the bottom boundary layer (BBL) parameterisation, first introduced by \citet{beckmann.doscher_JPO97}, 
     958A similar problem occurs in the $s$-coordinate when 
     959the thickness of the bottom level varies rapidly downstream of a sill 
     960\citep{willebrand.barnier.ea_PO01}, and the thickness of the plume is not resolved. 
     961 
     962The idea of the bottom boundary layer (BBL) parameterisation, first introduced by 
     963\citet{beckmann.doscher_JPO97}, 
    876964is to allow a direct communication between two adjacent bottom cells at different levels, 
    877965whenever the densest water is located above the less dense water. 
    878 The communication can be by a diffusive flux (diffusive BBL), an advective flux (advective BBL), or both. 
     966The communication can be by a diffusive flux (diffusive BBL), 
     967an advective flux (advective BBL), or both. 
    879968In the current implementation of the BBL, only the tracers are modified, not the velocities. 
    880 Furthermore, it only connects ocean bottom cells, and therefore does not include all the improvements introduced by 
    881 \citet{campin.goosse_T99}. 
     969Furthermore, it only connects ocean bottom cells, 
     970and therefore does not include all the improvements introduced by \citet{campin.goosse_T99}. 
    882971 
    883972%% ================================================================================================= 
     
    885974\label{subsec:TRA_bbl_diff} 
    886975 
    887 When applying sigma-diffusion (\np[=.true.]{ln_trabbl}{ln\_trabbl} and \np{nn_bbl_ldf}{nn\_bbl\_ldf} set to 1), 
     976When applying sigma-diffusion 
     977(\np[=.true.]{ln_trabbl}{ln\_trabbl} and \np{nn_bbl_ldf}{nn\_bbl\_ldf} set to 1), 
    888978the diffusive flux between two adjacent cells at the ocean floor is given by 
    889979\[ 
     
    891981  \vect F_\sigma = A_l^\sigma \, \nabla_\sigma T 
    892982\] 
    893 with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, and 
    894 $A_l^\sigma$ the lateral diffusivity in the BBL. 
     983with $\nabla_\sigma$ the lateral gradient operator taken between bottom cells, 
     984and $A_l^\sigma$ the lateral diffusivity in the BBL. 
    895985Following \citet{beckmann.doscher_JPO97}, the latter is prescribed with a spatial dependence, 
    896986\ie\ in the conditional form 
     
    900990      \begin{cases} 
    901991        A_{bbl} & \text{if~} \nabla_\sigma \rho \cdot \nabla H < 0 \\ 
    902         \\ 
    903         0      & \text{otherwise} \\ 
     992        0      & \text{otherwise} 
    904993      \end{cases} 
    905994\end{equation} 
    906 where $A_{bbl}$ is the BBL diffusivity coefficient, given by the namelist parameter \np{rn_ahtbbl}{rn\_ahtbbl} and 
     995where $A_{bbl}$ is the BBL diffusivity coefficient, 
     996given by the namelist parameter \np{rn_ahtbbl}{rn\_ahtbbl} and 
    907997usually set to a value much larger than the one used for lateral mixing in the open ocean. 
    908998The constraint in \autoref{eq:TRA_bbl_coef} implies that sigma-like diffusion only occurs when 
     
    9151005  \nabla_\sigma \rho / \rho = \alpha \, \nabla_\sigma T + \beta \, \nabla_\sigma S 
    9161006\] 
    917 where $\rho$, $\alpha$ and $\beta$ are functions of $\overline T^\sigma$, $\overline S^\sigma$ and 
    918 $\overline H^\sigma$, the along bottom mean temperature, salinity and depth, respectively. 
     1007where $\rho$, $\alpha$ and $\beta$ are functions of 
     1008$\overline T^\sigma$, $\overline S^\sigma$ and $\overline H^\sigma$, 
     1009the along bottom mean temperature, salinity and depth, respectively. 
    9191010 
    9201011%% ================================================================================================= 
     
    9271018%} 
    9281019 
    929 \begin{figure}[!t] 
     1020\begin{figure} 
    9301021  \centering 
    931   \includegraphics[width=0.66\textwidth]{Fig_BBL_adv} 
     1022  \includegraphics[width=0.33\textwidth]{TRA_BBL_adv} 
    9321023  \caption[Advective/diffusive bottom boundary layer]{ 
    9331024    Advective/diffusive Bottom Boundary Layer. 
     
    9461037%!!        i.e. transport proportional to the along-slope density gradient 
    9471038 
    948 %%%gmcomment   :  this section has to be really written 
    949  
    950 When applying an advective BBL (\np[=1..2]{nn_bbl_adv}{nn\_bbl\_adv}), an overturning circulation is added which 
    951 connects two adjacent bottom grid-points only if dense water overlies less dense water on the slope. 
     1039\cmtgm{This section has to be really written} 
     1040 
     1041When applying an advective BBL (\np[=1..2]{nn_bbl_adv}{nn\_bbl\_adv}), 
     1042an overturning circulation is added which connects two adjacent bottom grid-points only if 
     1043dense water overlies less dense water on the slope. 
    9521044The density difference causes dense water to move down the slope. 
    9531045 
    954 \np[=1]{nn_bbl_adv}{nn\_bbl\_adv}: 
    955 the downslope velocity is chosen to be the Eulerian ocean velocity just above the topographic step 
    956 (see black arrow in \autoref{fig:TRA_bbl}) \citep{beckmann.doscher_JPO97}. 
    957 It is a \textit{conditional advection}, that is, advection is allowed only 
    958 if dense water overlies less dense water on the slope (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$) and 
    959 if the velocity is directed towards greater depth (\ie\ $\vect U \cdot \nabla H > 0$). 
    960  
    961 \np[=2]{nn_bbl_adv}{nn\_bbl\_adv}: 
    962 the downslope velocity is chosen to be proportional to $\Delta \rho$, 
    963 the density difference between the higher cell and lower cell densities \citep{campin.goosse_T99}. 
    964 The advection is allowed only  if dense water overlies less dense water on the slope 
    965 (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$). 
    966 For example, the resulting transport of the downslope flow, here in the $i$-direction (\autoref{fig:TRA_bbl}), 
    967 is simply given by the following expression: 
    968 \[ 
    969   % \label{eq:TRA_bbl_Utr} 
    970   u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn}) 
    971 \] 
    972 where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as \np{rn_gambbl}{rn\_gambbl}, 
    973 a namelist parameter, and \textit{kup} and \textit{kdwn} are the vertical index of the higher and lower cells, 
    974 respectively. 
    975 The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity, 
    976 and because no direct estimation of this parameter is available, a uniform value has been assumed. 
    977 The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}. 
    978  
    979 Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using the upwind scheme. 
    980 Such a diffusive advective scheme has been chosen to mimic the entrainment between the downslope plume and 
    981 the surrounding water at intermediate depths. 
     1046\begin{description} 
     1047\item [{\np[=1]{nn_bbl_adv}{nn\_bbl\_adv}}] the downslope velocity is chosen to 
     1048  be the Eulerian ocean velocity just above the topographic step 
     1049  (see black arrow in \autoref{fig:TRA_bbl}) \citep{beckmann.doscher_JPO97}. 
     1050  It is a \textit{conditional advection}, that is, 
     1051  advection is allowed only if dense water overlies less dense water on the slope 
     1052  (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$) and if the velocity is directed towards greater depth 
     1053  (\ie\ $\vect U \cdot \nabla H > 0$). 
     1054\item [{\np[=2]{nn_bbl_adv}{nn\_bbl\_adv}}] the downslope velocity is chosen to be proportional to 
     1055  $\Delta \rho$, the density difference between the higher cell and lower cell densities 
     1056  \citep{campin.goosse_T99}. 
     1057  The advection is allowed only  if dense water overlies less dense water on the slope 
     1058  (\ie\ $\nabla_\sigma \rho \cdot \nabla H < 0$). 
     1059  For example, the resulting transport of the downslope flow, here in the $i$-direction 
     1060  (\autoref{fig:TRA_bbl}), is simply given by the following expression: 
     1061  \[ 
     1062    % \label{eq:TRA_bbl_Utr} 
     1063    u^{tr}_{bbl} = \gamma g \frac{\Delta \rho}{\rho_o} e_{1u} \, min ({e_{3u}}_{kup},{e_{3u}}_{kdwn}) 
     1064  \] 
     1065  where $\gamma$, expressed in seconds, is the coefficient of proportionality provided as 
     1066  \np{rn_gambbl}{rn\_gambbl}, a namelist parameter, and 
     1067  \textit{kup} and \textit{kdwn} are the vertical index of the higher and lower cells, respectively. 
     1068  The parameter $\gamma$ should take a different value for each bathymetric step, but for simplicity, 
     1069  and because no direct estimation of this parameter is available, a uniform value has been assumed. 
     1070  The possible values for $\gamma$ range between 1 and $10~s$ \citep{campin.goosse_T99}. 
     1071\end{description} 
     1072 
     1073Scalar properties are advected by this additional transport $(u^{tr}_{bbl},v^{tr}_{bbl})$ using 
     1074the upwind scheme. 
     1075Such a diffusive advective scheme has been chosen to mimic the entrainment between 
     1076the downslope plume and the surrounding water at intermediate depths. 
    9821077The entrainment is replaced by the vertical mixing implicit in the advection scheme. 
    9831078Let us consider as an example the case displayed in \autoref{fig:TRA_bbl} where 
    9841079the density at level $(i,kup)$ is larger than the one at level $(i,kdwn)$. 
    985 The advective BBL scheme modifies the tracer time tendency of the ocean cells near the topographic step by 
    986 the downslope flow \autoref{eq:TRA_bbl_dw}, the horizontal \autoref{eq:TRA_bbl_hor} and 
    987 the upward \autoref{eq:TRA_bbl_up} return flows as follows: 
    988 \begin{alignat}{3} 
     1080The advective BBL scheme modifies the tracer time tendency of 
     1081the ocean cells near the topographic step by the downslope flow \autoref{eq:TRA_bbl_dw}, 
     1082the horizontal \autoref{eq:TRA_bbl_hor} and the upward \autoref{eq:TRA_bbl_up} return flows as follows: 
     1083\begin{alignat}{5} 
    9891084  \label{eq:TRA_bbl_dw} 
    990   \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} 
    991                                 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\ 
     1085  \partial_t T^{do}_{kdw} &\equiv \partial_t T^{do}_{kdw} &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{kdw}} &&\lt( T^{sh}_{kup} - T^{do}_{kdw} \rt) \\ 
    9921086  \label{eq:TRA_bbl_hor} 
    993   \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} 
    994                                 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\ 
    995   % 
    996   \intertext{and for $k =kdw-1,\;..., \; kup$ :} 
    997   % 
     1087  \partial_t T^{sh}_{kup} &\equiv \partial_t T^{sh}_{kup} &&+ \frac{u^{tr}_{bbl}}{{b_t}^{sh}_{kup}} &&\lt( T^{do}_{kup} - T^{sh}_{kup} \rt) \\ 
     1088  \shortintertext{and for $k =kdw-1,\;..., \; kup$ :} 
    9981089  \label{eq:TRA_bbl_up} 
    999   \partial_t T^{do}_{k} &\equiv \partial_t S^{do}_{k} 
    1000                                 &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   &&\lt( T^{do}_{k +1} - T^{sh}_{k}   \rt) 
     1090  \partial_t T^{do}_{k}   &\equiv \partial_t S^{do}_{k}   &&+ \frac{u^{tr}_{bbl}}{{b_t}^{do}_{k}}   &&\lt( T^{do}_{k +1} - T^{sh}_{k}   \rt) 
    10011091\end{alignat} 
    10021092where $b_t$ is the $T$-cell volume. 
     
    10151105\end{listing} 
    10161106 
    1017 In some applications it can be useful to add a Newtonian damping term into the temperature and salinity equations: 
     1107In some applications it can be useful to add a Newtonian damping term into 
     1108the temperature and salinity equations: 
    10181109\begin{equation} 
    10191110  \label{eq:TRA_dmp} 
    1020   \begin{gathered} 
    1021     \pd[T]{t} = \cdots - \gamma (T - T_o) \\ 
    1022     \pd[S]{t} = \cdots - \gamma (S - S_o) 
    1023   \end{gathered} 
    1024 \end{equation} 
    1025 where $\gamma$ is the inverse of a time scale, and $T_o$ and $S_o$ are given temperature and salinity fields 
    1026 (usually a climatology). 
    1027 Options are defined through the  \nam{tra_dmp}{tra\_dmp} namelist variables. 
     1111    \pd[T]{t} = \cdots - \gamma (T - T_o) \qquad \pd[S]{t} = \cdots - \gamma (S - S_o) 
     1112\end{equation} 
     1113where $\gamma$ is the inverse of a time scale, 
     1114and $T_o$ and $S_o$ are given temperature and salinity fields (usually a climatology). 
     1115Options are defined through the \nam{tra_dmp}{tra\_dmp} namelist variables. 
    10281116The restoring term is added when the namelist parameter \np{ln_tradmp}{ln\_tradmp} is set to true. 
    1029 It also requires that both \np{ln_tsd_init}{ln\_tsd\_init} and \np{ln_tsd_dmp}{ln\_tsd\_dmp} are set to true in 
    1030 \nam{tsd}{tsd} namelist as well as \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures are correctly set 
     1117It also requires that both \np{ln_tsd_init}{ln\_tsd\_init} and 
     1118\np{ln_tsd_dmp}{ln\_tsd\_dmp} are set to true in \nam{tsd}{tsd} namelist as well as 
     1119\np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures are correctly set 
    10311120(\ie\ that $T_o$ and $S_o$ are provided in input files and read using \mdl{fldread}, 
    10321121see \autoref{subsec:SBC_fldread}). 
    1033 The restoring coefficient $\gamma$ is a three-dimensional array read in during the \rou{tra\_dmp\_init} routine. 
     1122The restoring coefficient $\gamma$ is a three-dimensional array read in during 
     1123the \rou{tra\_dmp\_init} routine. 
    10341124The file name is specified by the namelist variable \np{cn_resto}{cn\_resto}. 
    1035 The DMP\_TOOLS tool is provided to allow users to generate the netcdf file. 
     1125The \texttt{DMP\_TOOLS} are provided to allow users to generate the netcdf file. 
    10361126 
    10371127The two main cases in which \autoref{eq:TRA_dmp} is used are 
    1038 \textit{(a)} the specification of the boundary conditions along artificial walls of a limited domain basin and 
    1039 \textit{(b)} the computation of the velocity field associated with a given $T$-$S$ field 
    1040 (for example to build the initial state of a prognostic simulation, 
    1041 or to use the resulting velocity field for a passive tracer study). 
     1128\begin{enumerate*}[label=(\textit{\alph*})] 
     1129\item the specification of the boundary conditions along 
     1130  artificial walls of a limited domain basin and 
     1131\item the computation of the velocity field associated with a given $T$-$S$ field 
     1132  (for example to build the initial state of a prognostic simulation, 
     1133  or to use the resulting velocity field for a passive tracer study). 
     1134\end{enumerate*} 
    10421135The first case applies to regional models that have artificial walls instead of open boundaries. 
    1043 In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) whereas 
    1044 it is zero in the interior of the model domain. 
     1136In the vicinity of these walls, $\gamma$ takes large values (equivalent to a time scale of a few days) 
     1137whereas it is zero in the interior of the model domain. 
    10451138The second case corresponds to the use of the robust diagnostic method \citep{sarmiento.bryan_JGR82}. 
    10461139It allows us to find the velocity field consistent with the model dynamics whilst 
    10471140having a $T$, $S$ field close to a given climatological field ($T_o$, $S_o$). 
    10481141 
    1049 The robust diagnostic method is very efficient in preventing temperature drift in intermediate waters but 
    1050 it produces artificial sources of heat and salt within the ocean. 
     1142The robust diagnostic method is very efficient in preventing temperature drift in 
     1143intermediate waters but it produces artificial sources of heat and salt within the ocean. 
    10511144It also has undesirable effects on the ocean convection. 
    1052 It tends to prevent deep convection and subsequent deep-water formation, by stabilising the water column too much. 
    1053  
    1054 The namelist parameter \np{nn_zdmp}{nn\_zdmp} sets whether the damping should be applied in the whole water column or 
    1055 only below the mixed layer (defined either on a density or $S_o$ criterion). 
     1145It tends to prevent deep convection and subsequent deep-water formation, 
     1146by stabilising the water column too much. 
     1147 
     1148The namelist parameter \np{nn_zdmp}{nn\_zdmp} sets whether the damping should be applied in 
     1149the whole water column or only below the mixed layer (defined either on a density or $S_o$ criterion). 
    10561150It is common to set the damping to zero in the mixed layer as the adjustment time scale is short here 
    10571151\citep{madec.delecluse.ea_JPO96}. 
    10581152 
    1059 For generating \ifile{resto}, see the documentation for the DMP tool provided with the source code under 
    1060 \path{./tools/DMP_TOOLS}. 
     1153For generating \ifile{resto}, 
     1154see the documentation for the DMP tools provided with the source code under \path{./tools/DMP_TOOLS}. 
    10611155 
    10621156%% ================================================================================================= 
     
    10651159 
    10661160Options are defined through the \nam{dom}{dom} namelist variables. 
    1067 The general framework for tracer time stepping is a modified leap-frog scheme \citep{leclair.madec_OM09}, 
    1068 \ie\ a three level centred time scheme associated with a Asselin time filter (cf. \autoref{sec:TD_mLF}): 
     1161The general framework for tracer time stepping is a modified leap-frog scheme 
     1162\citep{leclair.madec_OM09}, \ie\ a three level centred time scheme associated with 
     1163a Asselin time filter (cf. \autoref{sec:TD_mLF}): 
    10691164\begin{equation} 
    10701165  \label{eq:TRA_nxt} 
    1071   \begin{alignedat}{3} 
     1166  \begin{alignedat}{5} 
    10721167    &(e_{3t}T)^{t + \rdt} &&= (e_{3t}T)_f^{t - \rdt} &&+ 2 \, \rdt \,e_{3t}^t \ \text{RHS}^t \\ 
    10731168    &(e_{3t}T)_f^t        &&= (e_{3t}T)^t            &&+ \, \gamma \, \lt[ (e_{3t}T)_f^{t - \rdt} - 2(e_{3t}T)^t + (e_{3t}T)^{t + \rdt} \rt] \\ 
     
    10751170  \end{alignedat} 
    10761171\end{equation} 
    1077 where RHS is the right hand side of the temperature equation, the subscript $f$ denotes filtered values, 
    1078 $\gamma$ is the Asselin coefficient, and $S$ is the total forcing applied on $T$ 
    1079 (\ie\ fluxes plus content in mass exchanges). 
    1080 $\gamma$ is initialized as \np{rn_atfp}{rn\_atfp} (\textbf{namelist} parameter). 
    1081 Its default value is \np[=10.e-3]{rn_atfp}{rn\_atfp}. 
     1172where RHS is the right hand side of the temperature equation, 
     1173the subscript $f$ denotes filtered values, $\gamma$ is the Asselin coefficient, 
     1174and $S$ is the total forcing applied on $T$ (\ie\ fluxes plus content in mass exchanges). 
     1175$\gamma$ is initialized as \np{rn_atfp}{rn\_atfp}, its default value is \forcode{10.e-3}. 
    10821176Note that the forcing correction term in the filter is not applied in linear free surface 
    10831177(\jp{ln\_linssh}\forcode{=.true.}) (see \autoref{subsec:TRA_sbc}). 
    1084 Not also that in constant volume case, the time stepping is performed on $T$, not on its content, $e_{3t}T$. 
    1085  
    1086 When the vertical mixing is solved implicitly, the update of the \textit{next} tracer fields is done in 
    1087 \mdl{trazdf} module. 
     1178Not also that in constant volume case, the time stepping is performed on $T$, 
     1179not on its content, $e_{3t}T$. 
     1180 
     1181When the vertical mixing is solved implicitly, 
     1182the update of the \textit{next} tracer fields is done in \mdl{trazdf} module. 
    10881183In this case only the swapping of arrays and the Asselin filtering is done in the \mdl{tranxt} module. 
    10891184 
    1090 In order to prepare for the computation of the \textit{next} time step, a swap of tracer arrays is performed: 
    1091 $T^{t - \rdt} = T^t$ and $T^t = T_f$. 
     1185In order to prepare for the computation of the \textit{next} time step, 
     1186a swap of tracer arrays is performed: $T^{t - \rdt} = T^t$ and $T^t = T_f$. 
    10921187 
    10931188%% ================================================================================================= 
     
    11051200\label{subsec:TRA_eos} 
    11061201 
    1107 The Equation Of Seawater (EOS) is an empirical nonlinear thermodynamic relationship linking seawater density, 
    1108 $\rho$, to a number of state variables, most typically temperature, salinity and pressure. 
     1202The \textbf{E}quation \textbf{O}f \textbf{S}eawater (EOS) is 
     1203an empirical nonlinear thermodynamic relationship linking 
     1204seawater density, $\rho$, to a number of state variables, 
     1205most typically temperature, salinity and pressure. 
    11091206Because density gradients control the pressure gradient force through the hydrostatic balance, 
    1110 the equation of state provides a fundamental bridge between the distribution of active tracers and 
    1111 the fluid dynamics. 
     1207the equation of state provides a fundamental bridge between 
     1208the distribution of active tracers and the fluid dynamics. 
    11121209Nonlinearities of the EOS are of major importance, in particular influencing the circulation through 
    11131210determination of the static stability below the mixed layer, 
    1114 thus controlling rates of exchange between the atmosphere and the ocean interior \citep{roquet.madec.ea_JPO15}. 
    1115 Therefore an accurate EOS based on either the 1980 equation of state (EOS-80, \cite{fofonoff.millard_bk83}) or 
    1116 TEOS-10 \citep{ioc.iapso_bk10} standards should be used anytime a simulation of the real ocean circulation is attempted 
     1211thus controlling rates of exchange between the atmosphere and the ocean interior 
    11171212\citep{roquet.madec.ea_JPO15}. 
     1213Therefore an accurate EOS based on either the 1980 equation of state 
     1214(EOS-80, \cite{fofonoff.millard_bk83}) or TEOS-10 \citep{ioc.iapso_bk10} standards should 
     1215be used anytime a simulation of the real ocean circulation is attempted \citep{roquet.madec.ea_JPO15}. 
    11181216The use of TEOS-10 is highly recommended because 
    1119 \textit{(i)}   it is the new official EOS, 
    1120 \textit{(ii)}  it is more accurate, being based on an updated database of laboratory measurements, and 
    1121 \textit{(iii)} it uses Conservative Temperature and Absolute Salinity (instead of potential temperature and 
    1122 practical salinity for EOS-80, both variables being more suitable for use as model variables 
    1123 \citep{ioc.iapso_bk10, graham.mcdougall_JPO13}. 
     1217\begin{enumerate*}[label=(\textit{\roman*})] 
     1218\item it is the new official EOS, 
     1219\item it is more accurate, being based on an updated database of laboratory measurements, and 
     1220\item it uses Conservative Temperature and Absolute Salinity 
     1221  (instead of potential temperature and practical salinity for EOS-80), 
     1222  both variables being more suitable for use as model variables 
     1223  \citep{ioc.iapso_bk10, graham.mcdougall_JPO13}. 
     1224\end{enumerate*} 
    11241225EOS-80 is an obsolescent feature of the \NEMO\ system, kept only for backward compatibility. 
    11251226For process studies, it is often convenient to use an approximation of the EOS. 
    11261227To that purposed, a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is also available. 
    11271228 
    1128 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, with $\rho_o$ a reference density. 
    1129 Called \textit{rau0} in the code, $\rho_o$ is set in \mdl{phycst} to a value of $1,026~Kg/m^3$. 
    1130 This is a sensible choice for the reference density used in a Boussinesq ocean climate model, as, 
    1131 with the exception of only a small percentage of the ocean, 
    1132 density in the World Ocean varies by no more than 2$\%$ from that value \citep{gill_bk82}. 
     1229In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed, 
     1230with $\rho_o$ a reference density. 
     1231Called \textit{rau0} in the code, 
     1232$\rho_o$ is set in \mdl{phycst} to a value of \texttt{1,026} $Kg/m^3$. 
     1233This is a sensible choice for the reference density used in a Boussinesq ocean climate model, 
     1234as, with the exception of only a small percentage of the ocean, 
     1235density in the World Ocean varies by no more than 2\% from that value \citep{gill_bk82}. 
    11331236 
    11341237Options which control the EOS used are defined through the \nam{eos}{eos} namelist variables. 
    11351238 
    11361239\begin{description} 
    1137 \item [{\np[=.true.]{ln_teos10}{ln\_teos10}}] the polyTEOS10-bsq equation of seawater \citep{roquet.madec.ea_OM15} is used. 
     1240\item [{\np[=.true.]{ln_teos10}{ln\_teos10}}] the polyTEOS10-bsq equation of seawater 
     1241  \citep{roquet.madec.ea_OM15} is used. 
    11381242  The accuracy of this approximation is comparable to the TEOS-10 rational function approximation, 
    1139   but it is optimized for a boussinesq fluid and the polynomial expressions have simpler and 
    1140   more computationally efficient expressions for their derived quantities which make them more adapted for 
    1141   use in ocean models. 
    1142   Note that a slightly higher precision polynomial form is now used replacement of 
    1143   the TEOS-10 rational function approximation for hydrographic data analysis \citep{ioc.iapso_bk10}. 
     1243  but it is optimized for a Boussinesq fluid and 
     1244  the polynomial expressions have simpler and more computationally efficient expressions for 
     1245  their derived quantities which make them more adapted for use in ocean models. 
     1246  Note that a slightly higher precision polynomial form is now used 
     1247  replacement of the TEOS-10 rational function approximation for hydrographic data analysis 
     1248  \citep{ioc.iapso_bk10}. 
    11441249  A key point is that conservative state variables are used: 
    1145   Absolute Salinity (unit: g/kg, notation: $S_A$) and Conservative Temperature (unit: \deg{C}, notation: $\Theta$). 
     1250  Absolute Salinity (unit: $g/kg$, notation: $S_A$) and 
     1251  Conservative Temperature (unit: $\deg{C}$, notation: $\Theta$). 
    11461252  The pressure in decibars is approximated by the depth in meters. 
    11471253  With TEOS10, the specific heat capacity of sea water, $C_p$, is a constant. 
    1148   It is set to $C_p = 3991.86795711963~J\,Kg^{-1}\,^{\circ}K^{-1}$, according to \citet{ioc.iapso_bk10}. 
     1254  It is set to $C_p$ = 3991.86795711963 $J.Kg^{-1}.\deg{K}^{-1}$, 
     1255  according to \citet{ioc.iapso_bk10}. 
    11491256  Choosing polyTEOS10-bsq implies that the state variables used by the model are $\Theta$ and $S_A$. 
    1150   In particular, the initial state deined by the user have to be given as \textit{Conservative} Temperature and 
    1151   \textit{Absolute} Salinity. 
     1257  In particular, the initial state defined by the user have to be given as 
     1258  \textit{Conservative} Temperature and \textit{Absolute} Salinity. 
    11521259  In addition, when using TEOS10, the Conservative SST is converted to potential SST prior to 
    11531260  either computing the air-sea and ice-sea fluxes (forced mode) or 
    11541261  sending the SST field to the atmosphere (coupled mode). 
    11551262\item [{\np[=.true.]{ln_eos80}{ln\_eos80}}] the polyEOS80-bsq equation of seawater is used. 
    1156   It takes the same polynomial form as the polyTEOS10, but the coefficients have been optimized to 
    1157   accurately fit EOS80 (Roquet, personal comm.). 
     1263  It takes the same polynomial form as the polyTEOS10, 
     1264  but the coefficients have been optimized to accurately fit EOS80 (Roquet, personal comm.). 
    11581265  The state variables used in both the EOS80 and the ocean model are: 
    1159   the Practical Salinity ((unit: psu, notation: $S_p$)) and 
    1160   Potential Temperature (unit: $^{\circ}C$, notation: $\theta$). 
     1266  the Practical Salinity (unit: $psu$, notation: $S_p$) and 
     1267  Potential Temperature (unit: $\deg{C}$, notation: $\theta$). 
    11611268  The pressure in decibars is approximated by the depth in meters. 
    1162   With thsi EOS, the specific heat capacity of sea water, $C_p$, is a function of temperature, salinity and 
    1163   pressure \citep{fofonoff.millard_bk83}. 
     1269  With EOS, the specific heat capacity of sea water, $C_p$, is a function of 
     1270  temperature, salinity and pressure \citep{fofonoff.millard_bk83}. 
    11641271  Nevertheless, a severe assumption is made in order to have a heat content ($C_p T_p$) which 
    11651272  is conserved by the model: $C_p$ is set to a constant value, the TEOS10 value. 
    1166 \item [{\np[=.true.]{ln_seos}{ln\_seos}}] a simplified EOS (S-EOS) inspired by \citet{vallis_bk06} is chosen, 
    1167   the coefficients of which has been optimized to fit the behavior of TEOS10 
    1168   (Roquet, personal comm.) (see also \citet{roquet.madec.ea_JPO15}). 
     1273\item [{\np[=.true.]{ln_seos}{ln\_seos}}] a simplified EOS (S-EOS) inspired by 
     1274  \citet{vallis_bk06} is chosen, 
     1275  the coefficients of which has been optimized to fit the behavior of TEOS10 (Roquet, personal comm.) 
     1276  (see also \citet{roquet.madec.ea_JPO15}). 
    11691277  It provides a simplistic linear representation of both cabbeling and thermobaricity effects which 
    11701278  is enough for a proper treatment of the EOS in theoretical studies \citep{roquet.madec.ea_JPO15}. 
    1171   With such an equation of state there is no longer a distinction between 
    1172   \textit{conservative} and \textit{potential} temperature, 
    1173   as well as between \textit{absolute} and \textit{practical} salinity. 
     1279  With such an equation of state there is no longer a distinction between \textit{conservative} and 
     1280  \textit{potential} temperature, as well as between \textit{absolute} and 
     1281  \textit{practical} salinity. 
    11741282  S-EOS takes the following expression: 
    1175  
    11761283  \begin{gather*} 
    11771284    % \label{eq:TRA_S-EOS} 
    1178     \begin{alignedat}{2} 
    1179     &d_a(T,S,z) = \frac{1}{\rho_o} \big[ &- a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * &T_a \big. \\ 
    1180     &                                    &+ b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * &S_a       \\ 
    1181     &                              \big. &- \nu \;                           T_a                  &S_a \big] \\ 
    1182     \end{alignedat} 
    1183     \\ 
     1285    d_a(T,S,z) = \frac{1}{\rho_o} \big[ - a_0 \; ( 1 + 0.5 \; \lambda_1 \; T_a + \mu_1 \; z ) * T_a \big. 
     1286                                        + b_0 \; ( 1 - 0.5 \; \lambda_2 \; S_a - \mu_2 \; z ) * S_a 
     1287                                  \big. - \nu \;                           T_a                  S_a \big] \\ 
    11841288    \text{with~} T_a = T - 10 \, ; \, S_a = S - 35 \, ; \, \rho_o = 1026~Kg/m^3 
    11851289  \end{gather*} 
    1186   where the computer name of the coefficients as well as their standard value are given in \autoref{tab:TRA_SEOS}. 
     1290  where the computer name of the coefficients as well as their standard value are given in 
     1291  \autoref{tab:TRA_SEOS}. 
    11871292  In fact, when choosing S-EOS, various approximation of EOS can be specified simply by 
    11881293  changing the associated coefficients. 
    1189   Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$ remove thermobaric effect from S-EOS. 
    1190   setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$ remove cabbeling effect from 
    1191   S-EOS. 
     1294  Setting to zero the two thermobaric coefficients $(\mu_1,\mu_2)$ 
     1295  remove thermobaric effect from S-EOS. 
     1296  Setting to zero the three cabbeling coefficients $(\lambda_1,\lambda_2,\nu)$ 
     1297  remove   cabbeling effect from S-EOS. 
    11921298  Keeping non-zero value to $a_0$ and $b_0$ provide a linear EOS function of T and S. 
    11931299\end{description} 
    11941300 
    1195 \begin{table}[!tb] 
     1301\begin{table} 
    11961302  \centering 
    11971303  \begin{tabular}{|l|l|l|l|} 
    11981304    \hline 
    1199     coeff.     & computer name   & S-EOS           & description                      \\ 
     1305    coeff.      & computer name                & S-EOS            & description                     \\ 
    12001306    \hline 
    1201     $a_0$       & \np{rn_a0}{rn\_a0}     & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 
     1307    $a_0      $ & \np{rn_a0}{rn\_a0}           & $1.6550~10^{-1}$ & linear thermal expansion coeff. \\ 
    12021308    \hline 
    1203     $b_0$         & \np{rn_b0}{rn\_b0}       & $7.6554~10^{-1}$ & linear haline  expansion coeff. \\ 
     1309    $b_0      $ & \np{rn_b0}{rn\_b0}           & $7.6554~10^{-1}$ & linear haline  expansion coeff. \\ 
    12041310    \hline 
    1205     $\lambda_1$   & \np{rn_lambda1}{rn\_lambda1}& $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$       \\ 
     1311    $\lambda_1$ & \np{rn_lambda1}{rn\_lambda1} & $5.9520~10^{-2}$ & cabbeling coeff. in $T^2$       \\ 
    12061312    \hline 
    1207     $\lambda_2$   & \np{rn_lambda2}{rn\_lambda2}& $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$       \\ 
     1313    $\lambda_2$ & \np{rn_lambda2}{rn\_lambda2} & $5.4914~10^{-4}$ & cabbeling coeff. in $S^2$       \\ 
    12081314    \hline 
    1209     $\nu$       & \np{rn_nu}{rn\_nu}     & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$     \\ 
     1315    $\nu      $ & \np{rn_nu}{rn\_nu}           & $2.4341~10^{-3}$ & cabbeling coeff. in $T \, S$    \\ 
    12101316    \hline 
    1211     $\mu_1$     & \np{rn_mu1}{rn\_mu1}    & $1.4970~10^{-4}$ & thermobaric coeff. in T         \\ 
     1317    $\mu_1    $ & \np{rn_mu1}{rn\_mu1}         & $1.4970~10^{-4}$ & thermobaric coeff. in T         \\ 
    12121318    \hline 
    1213     $\mu_2$     & \np{rn_mu2}{rn\_mu2}    & $1.1090~10^{-5}$ & thermobaric coeff. in S         \\ 
     1319    $\mu_2    $ & \np{rn_mu2}{rn\_mu2}         & $1.1090~10^{-5}$ & thermobaric coeff. in S         \\ 
    12141320    \hline 
    12151321  \end{tabular} 
     
    12221328\label{subsec:TRA_bn2} 
    12231329 
    1224 An accurate computation of the ocean stability (i.e. of $N$, the brunt-V\"{a}is\"{a}l\"{a} frequency) is of 
    1225 paramount importance as determine the ocean stratification and is used in several ocean parameterisations 
     1330An accurate computation of the ocean stability (i.e. of $N$, the Brunt-V\"{a}is\"{a}l\"{a} frequency) is of paramount importance as determine the ocean stratification and 
     1331is used in several ocean parameterisations 
    12261332(namely TKE, GLS, Richardson number dependent vertical diffusion, enhanced vertical diffusion, 
    12271333non-penetrative convection, tidal mixing  parameterisation, iso-neutral diffusion). 
     
    12351341where $(T,S) = (\Theta,S_A)$ for TEOS10, $(\theta,S_p)$ for TEOS-80, or $(T,S)$ for S-EOS, and, 
    12361342$\alpha$ and $\beta$ are the thermal and haline expansion coefficients. 
    1237 The coefficients are a polynomial function of temperature, salinity and depth which expression depends on 
    1238 the chosen EOS. 
     1343The coefficients are a polynomial function of temperature, salinity and depth which 
     1344expression depends on the chosen EOS. 
    12391345They are computed through \textit{eos\_rab}, a \fortran\ function that can be found in \mdl{eosbn2}. 
    12401346 
     
    12461352\begin{equation} 
    12471353  \label{eq:TRA_eos_fzp} 
    1248   \begin{split} 
    1249     &T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\ 
    1250     &\text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \\ 
    1251     &\text{and~} d = -7.53~10^{-3} 
    1252     \end{split} 
     1354  \begin{gathered} 
     1355    T_f (S,p) = \lt( a + b \, \sqrt{S} + c \, S \rt) \, S + d \, p \\ 
     1356    \text{where~} a = -0.0575, \, b = 1.710523~10^{-3}, \, c = -2.154996~10^{-4} \text{and~} d = -7.53~10^{-3} 
     1357    \end{gathered} 
    12531358\end{equation} 
    12541359 
     
    12691374\label{sec:TRA_zpshde} 
    12701375 
    1271 \gmcomment{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 
     1376\cmtgm{STEVEN: to be consistent with earlier discussion of differencing and averaging operators, 
    12721377I've changed "derivative" to "difference" and "mean" to "average"} 
    12731378 
    1274 With partial cells (\np[=.true.]{ln_zps}{ln\_zps}) at bottom and top (\np[=.true.]{ln_isfcav}{ln\_isfcav}), 
     1379With partial cells (\np[=.true.]{ln_zps}{ln\_zps}) at bottom and top 
     1380(\np[=.true.]{ln_isfcav}{ln\_isfcav}), 
    12751381in general, tracers in horizontally adjacent cells live at different depths. 
    1276 Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) and 
    1277 the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 
    1278 The partial cell properties at the top (\np[=.true.]{ln_isfcav}{ln\_isfcav}) are computed in the same way as 
    1279 for the bottom. 
     1382Horizontal gradients of tracers are needed for horizontal diffusion 
     1383(\mdl{traldf} module) and the hydrostatic pressure gradient calculations (\mdl{dynhpg} module). 
     1384The partial cell properties at the top (\np[=.true.]{ln_isfcav}{ln\_isfcav}) are computed in 
     1385the same way as for the bottom. 
    12801386So, only the bottom interpolation is explained below. 
    12811387 
     
    12831389a linear interpolation in the vertical is used to approximate the deeper tracer as if 
    12841390it actually lived at the depth of the shallower tracer point (\autoref{fig:TRA_Partial_step_scheme}). 
    1285 For example, for temperature in the $i$-direction the needed interpolated temperature, $\widetilde T$, is: 
    1286  
    1287 \begin{figure}[!p] 
     1391For example, for temperature in the $i$-direction the needed interpolated temperature, 
     1392$\widetilde T$, is: 
     1393 
     1394\begin{figure} 
    12881395  \centering 
    1289   \includegraphics[width=0.66\textwidth]{Fig_partial_step_scheme} 
     1396  \includegraphics[width=0.33\textwidth]{TRA_partial_step_scheme} 
    12901397  \caption[Discretisation of the horizontal difference and average of tracers in 
    12911398  the $z$-partial step coordinate]{ 
     
    12941401    the case $(e3w_k^{i + 1} - e3w_k^i) > 0$. 
    12951402    A linear interpolation is used to estimate $\widetilde T_k^{i + 1}$, 
    1296     the tracer value at the depth of the shallower tracer point of 
    1297     the two adjacent bottom $T$-points. 
     1403    the tracer value at the depth of the shallower tracer point of the two adjacent bottom $T$-points. 
    12981404    The horizontal difference is then given by: 
    1299     $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and 
    1300     the average by: 
     1405    $\delta_{i + 1/2} T_k = \widetilde T_k^{\, i + 1} -T_k^{\, i}$ and the average by: 
    13011406    $\overline T_k^{\, i + 1/2} = (\widetilde T_k^{\, i + 1/2} - T_k^{\, i}) / 2$.} 
    13021407  \label{fig:TRA_Partial_step_scheme} 
    13031408\end{figure} 
     1409 
    13041410\[ 
    13051411  \widetilde T = \lt\{ 
    13061412    \begin{alignedat}{2} 
    13071413      &T^{\, i + 1} &-\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt) }{ e_{3w}^{i + 1} } \; \delta_k T^{i + 1} 
    1308       & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\ \\ 
     1414      & \quad \text{if $e_{3w}^{i + 1} \geq e_{3w}^i$} \\ 
    13091415      &T^{\, i}     &+\frac{ \lt( e_{3w}^{i + 1} -e_{3w}^i \rt )}{e_{3w}^i       } \; \delta_k T^{i + 1} 
    13101416      & \quad \text{if $e_{3w}^{i + 1} <    e_{3w}^i$} 
     
    13121418  \rt. 
    13131419\] 
    1314 and the resulting forms for the horizontal difference and the horizontal average value of $T$ at a $U$-point are: 
     1420and the resulting forms for the horizontal difference and the horizontal average value of 
     1421$T$ at a $U$-point are: 
    13151422\begin{equation} 
    13161423  \label{eq:TRA_zps_hde} 
     
    13181425    \delta_{i + 1/2} T       &= 
    13191426    \begin{cases} 
    1320                                 \widetilde T - T^i          & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
    1321                                 \\ 
    1322                                 T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} <    e_{3w}^i 
    1323     \end{cases} 
    1324     \\ 
     1427      \widetilde T - T^i          & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
     1428      T^{\, i + 1} - \widetilde T & \text{if~} e_{3w}^{i + 1} <    e_{3w}^i 
     1429    \end{cases} \\ 
    13251430    \overline T^{\, i + 1/2} &= 
    13261431    \begin{cases} 
    1327                                 (\widetilde T - T^{\, i}   ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
    1328                                 \\ 
    1329                                 (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} <   e_{3w}^i 
     1432      (\widetilde T - T^{\, i}    ) / 2 & \text{if~} e_{3w}^{i + 1} \geq e_{3w}^i \\ 
     1433      (T^{\, i + 1} - \widetilde T) / 2 & \text{if~} e_{3w}^{i + 1} <   e_{3w}^i 
    13301434    \end{cases} 
    13311435  \end{split} 
     
    13341438The computation of horizontal derivative of tracers as well as of density is performed once for all at 
    13351439each time step in \mdl{zpshde} module and stored in shared arrays to be used when needed. 
    1336 It has to be emphasized that the procedure used to compute the interpolated density, $\widetilde \rho$, 
    1337 is not the same as that used for $T$ and $S$. 
    1338 Instead of forming a linear approximation of density, we compute $\widetilde \rho$ from the interpolated values of 
    1339 $T$ and $S$, and the pressure at a $u$-point 
     1440It has to be emphasized that the procedure used to compute the interpolated density, 
     1441$\widetilde \rho$, is not the same as that used for $T$ and $S$. 
     1442Instead of forming a linear approximation of density, 
     1443we compute $\widetilde \rho$ from the interpolated values of $T$ and $S$, 
     1444and the pressure at a $u$-point 
    13401445(in the equation of state pressure is approximated by depth, see \autoref{subsec:TRA_eos}): 
    13411446\[ 
     
    13451450 
    13461451This is a much better approximation as the variation of $\rho$ with depth (and thus pressure) 
    1347 is highly non-linear with a true equation of state and thus is badly approximated with a linear interpolation. 
    1348 This approximation is used to compute both the horizontal pressure gradient (\autoref{sec:DYN_hpg}) and 
    1349 the slopes of neutral surfaces (\autoref{sec:LDF_slp}). 
    1350  
    1351 Note that in almost all the advection schemes presented in this Chapter, 
     1452is highly non-linear with a true equation of state and thus is badly approximated with 
     1453a linear interpolation. 
     1454This approximation is used to compute both the horizontal pressure gradient (\autoref{sec:DYN_hpg}) 
     1455and the slopes of neutral surfaces (\autoref{sec:LDF_slp}). 
     1456 
     1457Note that in almost all the advection schemes presented in this chapter, 
    13521458both averaging and differencing operators appear. 
    13531459Yet \autoref{eq:TRA_zps_hde} has not been used in these schemes: 
     
    13561462The main motivation is to preserve the domain averaged mean variance of the advected field when 
    13571463using the $2^{nd}$ order centred scheme. 
    1358 Sensitivity of the advection schemes to the way horizontal averages are performed in the vicinity of 
    1359 partial cells should be further investigated in the near future. 
    1360 \gmcomment{gm :   this last remark has to be done} 
    1361  
    1362 \onlyinsubfile{\input{../../global/epilogue}} 
     1464Sensitivity of the advection schemes to the way horizontal averages are performed in 
     1465the vicinity of partial cells should be further investigated in the near future. 
     1466\cmtgm{gm :   this last remark has to be done} 
     1467 
     1468\subinc{\input{../../global/epilogue}} 
    13631469 
    13641470\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_ZDF.tex

    r11599 r11954  
    2828\clearpage 
    2929 
    30 %gm% Add here a small introduction to ZDF and naming of the different physics (similar to what have been written for TRA and DYN. 
     30\cmtgm{ Add here a small introduction to ZDF and naming of the different physics 
     31(similar to what have been written for TRA and DYN).} 
    3132 
    3233%% ================================================================================================= 
     
    248249\begin{figure}[!t] 
    249250  \centering 
    250   \includegraphics[width=0.66\textwidth]{Fig_mixing_length} 
     251  \includegraphics[width=0.66\textwidth]{ZDF_mixing_length} 
    251252  \caption[Mixing length computation]{Illustration of the mixing length computation} 
    252253  \label{fig:ZDF_mixing_length} 
     
    533534 in \citet{reffray.guillaume.ea_GMD15} for the \NEMO\ model. 
    534535 
    535 %% ================================================================================================= 
    536 \subsection[OSM: OSMosis boundary layer scheme (\forcode{ln_zdfosm})]{OSM: OSMosis boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})} 
     536% ------------------------------------------------------------------------------------------------------------- 
     537%        OSM OSMOSIS BL Scheme 
     538% ------------------------------------------------------------------------------------------------------------- 
     539\subsection[OSM: OSMOSIS boundary layer scheme (\forcode{ln_zdfosm = .true.})] 
     540{OSM: OSMOSIS boundary layer scheme (\protect\np{ln_zdfosm}{ln\_zdfosm})} 
    537541\label{subsec:ZDF_osm} 
    538542 
     
    543547\end{listing} 
    544548 
    545 The OSMOSIS turbulent closure scheme is based on......   TBC 
     549%-------------------------------------------------------------------------------------------------------------- 
     550\paragraph{Namelist choices} 
     551Most of the namelist options refer to how to specify the Stokes 
     552surface drift and penetration depth. There are three options: 
     553\begin{description} 
     554  \item \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave} Default value in \texttt{namelist\_ref}. In this case the Stokes drift is 
     555      assumed to be parallel to the surface wind stress, with 
     556      magnitude consistent with a constant turbulent Langmuir number 
     557    $\mathrm{La}_t=$ \protect\np{rn_m_la} {rn\_m\_la} i.e.\ 
     558    $u_{s0}=\tau/(\mathrm{La}_t^2\rho_0)$.  Default value of 
     559    \protect\np{rn_m_la}{rn\_m\_la} is 0.3. The Stokes penetration 
     560      depth $\delta = $ \protect\np{rn_osm_dstokes}{rn\_osm\_dstokes}; this has default value 
     561      of 5~m. 
     562 
     563  \item \protect\np[=1]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 
     564      assumed to be parallel to the surface wind stress, with 
     565      magnitude as in the classical Pierson-Moskowitz wind-sea 
     566      spectrum.  Significant wave height and 
     567      wave-mean period taken from this spectrum are used to calculate the Stokes penetration 
     568      depth, following the approach set out in  \citet{breivik.janssen.ea_JPO14}. 
     569 
     570    \item \protect\np[=2]{nn_osm_wave}{nn\_osm\_wave} In this case the Stokes drift is 
     571      taken from  ECMWF wave model output, though only the component parallel 
     572      to the wind stress is retained. Significant wave height and 
     573      wave-mean period from ECMWF wave model output are used to calculate the Stokes penetration 
     574      depth, again following \citet{breivik.janssen.ea_JPO14}. 
     575 
     576    \end{description} 
     577 
     578    Others refer to the treatment of diffusion and viscosity beneath 
     579    the surface boundary layer: 
     580\begin{description} 
     581   \item \protect\np{ln_kpprimix} {ln\_kpprimix}  Default is \np{.true.}. Switches on KPP-style Ri \#-dependent 
     582     mixing below the surface boundary layer. If this is set 
     583     \texttt{.true.}  the following variable settings are honoured: 
     584    \item \protect\np{rn_riinfty}{rn\_riinfty} Critical value of local Ri \# below which 
     585      shear instability increases vertical mixing from background value. 
     586    \item \protect\np{rn_difri} {rn\_difri} Maximum value of Ri \#-dependent mixing at $\mathrm{Ri}=0$. 
     587    \item \protect\np{ln_convmix}{ln\_convmix} If \texttt{.true.} then, where water column is unstable, specify 
     588       diffusivity equal to \protect\np{rn_dif_conv}{rn\_dif\_conv} (default value is 1 m~s$^{-2}$). 
     589 \end{description} 
     590 Diagnostic output is controlled by: 
     591  \begin{description} 
     592    \item \protect\np{ln_dia_osm}{ln\_dia\_osm} Default is \np{.false.}; allows XIOS output of OSMOSIS internal fields. 
     593  \end{description} 
     594Obsolete namelist parameters include: 
     595\begin{description} 
     596   \item \protect\np{ln_use_osm_la}\np{ln\_use\_osm\_la} With \protect\np[=0]{nn_osm_wave}{nn\_osm\_wave}, 
     597      \protect\np{rn_osm_dstokes} {rn\_osm\_dstokes} is always used to specify the Stokes 
     598      penetration depth. 
     599   \item \protect\np{nn_ave} {nn\_ave} Choice of averaging method for KPP-style Ri \# 
     600      mixing. Not taken account of. 
     601   \item \protect\np{rn_osm_hbl0} {rn\_osm\_hbl0} Depth of initial boundary layer is now set 
     602     by a density criterion similar to that used in calculating \emph{hmlp} (output as \texttt{mldr10\_1}) in \mdl{zdfmxl}. 
     603\end{description} 
     604 
     605\subsubsection{Summary} 
     606Much of the time the turbulent motions in the ocean surface boundary 
     607layer (OSBL) are not given by 
     608classical shear turbulence. Instead they are in a regime known as 
     609`Langmuir turbulence',  dominated by an 
     610interaction between the currents and the Stokes drift of the surface waves \citep[e.g.][]{mcwilliams.ea_JFM97}. 
     611This regime is characterised by strong vertical turbulent motion, and appears when the surface Stokes drift $u_{s0}$ is much greater than the friction velocity $u_{\ast}$. More specifically Langmuir turbulence is thought to be crucial where the turbulent Langmuir number $\mathrm{La}_{t}=(u_{\ast}/u_{s0}) > 0.4$. 
     612 
     613The OSMOSIS model is fundamentally based on results of Large Eddy 
     614Simulations (LES) of Langmuir turbulence and aims to fully describe 
     615this Langmuir regime. The description in this section is of necessity incomplete and further details are available in Grant. A (2019); in prep. 
     616 
     617The OSMOSIS turbulent closure scheme is a similarity-scale scheme in 
     618the same spirit as the K-profile 
     619parameterization (KPP) scheme of \citet{large.ea_RG97}. 
     620A specified shape of diffusivity, scaled by the (OSBL) depth 
     621$h_{\mathrm{BL}}$ and a turbulent velocity scale, is imposed throughout the 
     622boundary layer 
     623$-h_{\mathrm{BL}}<z<\eta$. The turbulent closure model 
     624also includes fluxes of tracers and momentum that are``non-local'' (independent of the local property gradient). 
     625 
     626Rather than the OSBL 
     627depth being diagnosed in terms of a bulk Richardson number criterion, 
     628as in KPP, it is set by a prognostic equation that is informed by 
     629energy budget considerations reminiscent of the classical mixed layer 
     630models of \citet{kraus.turner_tellus67}. 
     631The model also includes an explicit parametrization of the structure 
     632of the pycnocline (the stratified region at the bottom of the OSBL). 
     633 
     634Presently, mixing below the OSBL is handled by the Richardson 
     635number-dependent mixing scheme used in \citet{large.ea_RG97}. 
     636 
     637Convective parameterizations such as described in \ref{sec:ZDF_conv} 
     638below should not be used with the OSMOSIS-OBL model: instabilities 
     639within the OSBL are part of the model, while instabilities below the 
     640ML are handled by the Ri \# dependent scheme. 
     641 
     642\subsubsection{Depth and velocity scales} 
     643The model supposes a boundary layer of thickness $h_{\mathrm{bl}}$ enclosing a well-mixed layer of thickness $h_{\mathrm{ml}}$ and a relatively thin pycnocline at the base of thickness $\Delta h$; Fig.~\ref{fig: OSBL_structure} shows typical (a) buoyancy structure and (b) turbulent buoyancy flux profile for the unstable boundary layer (losing buoyancy at the surface; e.g.\ cooling). 
     644\begin{figure}[!t] 
     645  \begin{center} 
     646    %\includegraphics[width=0.7\textwidth]{ZDF_OSM_structure_of_OSBL} 
     647    \caption{ 
     648      \protect\label{fig: OSBL_structure} 
     649     The structure of the entraining  boundary layer. (a) Mean buoyancy profile. (b) Profile of the buoyancy flux. 
     650    } 
     651  \end{center} 
     652\end{figure} 
     653The pycnocline in the OSMOSIS scheme is assumed to have a finite thickness, and may include a number of model levels. This means that the OSMOSIS scheme must parametrize both the thickness of the pycnocline, and the turbulent fluxes within the pycnocline. 
     654 
     655Consideration of the power input by wind acting on the Stokes drift suggests that the Langmuir turbulence has velocity scale: 
     656\begin{equation}\label{eq:w_La} 
     657w_{*L}= \left(u_*^2 u_{s\,0}\right)^{1/3}; 
     658\end{equation} 
     659but at times the Stokes drift may be weak due to e.g.\ ice cover, short fetch, misalignment with the surface stress, etc.\ so  a composite velocity scale is assumed for the stable (warming) boundary layer: 
     660\begin{equation}\label{eq:composite-nu} 
     661  \nu_{\ast}= \left\{ u_*^3 \left[1-\exp(-.5 \mathrm{La}_t^2)\right]+w_{*L}^3\right\}^{1/3}. 
     662\end{equation} 
     663For the unstable boundary layer this is merged with the standard convective velocity scale $w_{*C}=\left(\overline{w^\prime b^\prime}_0 \,h_\mathrm{ml}\right)^{1/3}$, where $\overline{w^\prime b^\prime}_0$ is the upwards surface buoyancy flux, to give: 
     664\begin{equation}\label{eq:vel-scale-unstable} 
     665\omega_* = \left(\nu_*^3 + 0.5 w_{*C}^3\right)^{1/3}. 
     666\end{equation} 
     667 
     668\subsubsection{The flux gradient model} 
     669The flux-gradient relationships used in the OSMOSIS scheme take the form: 
     670% 
     671\begin{equation}\label{eq:flux-grad-gen} 
     672\overline{w^\prime\chi^\prime}=-K\frac{\partial\overline{\chi}}{\partial z} + N_{\chi,s} +N_{\chi,b} +N_{\chi,t}, 
     673\end{equation} 
     674% 
     675where $\chi$ is a general variable and $N_{\chi,s}, N_{\chi,b} \mathrm{and} N_{\chi,t}$  are the non-gradient terms, and represent the effects of the different terms in the turbulent flux-budget on the transport of $\chi$. $N_{\chi,s}$ represents the effects that the Stokes shear has on the transport of $\chi$, $N_{\chi,b}$  the effect of buoyancy, and $N_{\chi,t}$ the effect of the turbulent transport.  The same general form for the flux-gradient relationship is used to parametrize the transports of momentum, heat and salinity. 
     676 
     677In terms of the non-dimensionalized depth variables 
     678% 
     679\begin{equation}\label{eq:sigma} 
     680\sigma_{\mathrm{ml}}= -z/h_{\mathrm{ml}}; \;\sigma_{\mathrm{bl}}= -z/h_{\mathrm{bl}}, 
     681\end{equation} 
     682% 
     683in unstable conditions the eddy diffusivity ($K_d$) and eddy viscosity ($K_\nu$) profiles are parametrized as: 
     684% 
     685\begin{align}\label{eq:diff-unstable} 
     686K_d=&0.8\, \omega_*\, h_{\mathrm{ml}} \, \sigma_{\mathrm{ml}} \left(1-\beta_d \sigma_{\mathrm{ml}}\right)^{3/2} 
     687\\\label{eq:visc-unstable} 
     688K_\nu =& 0.3\, \omega_* \,h_{\mathrm{ml}}\, \sigma_{\mathrm{ml}} \left(1-\beta_\nu \sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right) 
     689\end{align} 
     690% 
     691where $\beta_d$ and $\beta_\nu$ are parameters that are determined by matching Eqs \ref{eq:diff-unstable} and \ref{eq:visc-unstable} to the eddy diffusivity and viscosity at the base of the well-mixed layer, given by 
     692% 
     693\begin{equation}\label{eq:diff-wml-base} 
     694K_{d,\mathrm{ml}}=K_{\nu,\mathrm{ml}}=\,0.16\,\omega_* \Delta h. 
     695\end{equation} 
     696% 
     697For stable conditions the eddy diffusivity/viscosity profiles are given by: 
     698% 
     699\begin{align}\label{diff-stable} 
     700K_d= & 0.75\,\, \nu_*\, h_{\mathrm{ml}}\,\,  \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right]\sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)^{3/2} \\\label{eq:visc-stable} 
     701K_\nu = & 0.375\,\,  \nu_*\, h_{\mathrm{ml}} \,\, \exp\left[-2.8 \left(h_{\mathrm{bl}}/L_L\right)^2\right] \sigma_{\mathrm{ml}} \left(1-\sigma_{\mathrm{ml}}\right)\left(1-\tfrac{1}{2}\sigma_{\mathrm{ml}}^2\right). 
     702\end{align} 
     703% 
     704The shape of the eddy viscosity and diffusivity profiles is the same as the shape in the unstable OSBL. The eddy diffusivity/viscosity depends on the stability parameter $h_{\mathrm{bl}}/{L_L}$ where $ L_L$ is analogous to the Obukhov length, but for Langmuir turbulence: 
     705\begin{equation}\label{eq:L_L} 
     706  L_L=-w_{*L}^3/\left<\overline{w^\prime b^\prime}\right>_L, 
     707\end{equation} 
     708with the mean turbulent buoyancy flux averaged over the boundary layer given in terms of its surface value $\overline{w^\prime b^\prime}_0$ and (downwards) )solar irradiance $I(z)$ by 
     709\begin{equation} \label{eq:stable-av-buoy-flux} 
     710\left<\overline{w^\prime b^\prime}\right>_L = \tfrac{1}{2} {\overline{w^\prime b^\prime}}_0-g\alpha_E\left[\tfrac{1}{2}(I(0)+I(-h))-\left<I\right>\right]. 
     711\end{equation} 
     712% 
     713In unstable conditions the eddy diffusivity and viscosity depend on stability through the velocity scale $\omega_*$, which depends on the two velocity scales $\nu_*$ and $w_{*C}$. 
     714 
     715Details of the non-gradient terms in \eqref{eq:flux-grad-gen} and of the fluxes within the pycnocline $-h_{\mathrm{bl}}<z<h_{\mathrm{ml}}$ can be found in Grant (2019). 
     716 
     717\subsubsection{Evolution of the boundary layer depth} 
     718 
     719The prognostic equation for the depth of the neutral/unstable boundary layer is given by \citep{grant+etal18}, 
     720 
     721\begin{equation} \label{eq:dhdt-unstable} 
     722%\frac{\partial h_\mathrm{bl}}{\partial t} + \mathbf{U}_b\cdot\nabla h_\mathrm{bl}= W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}} 
     723\frac{\partial h_\mathrm{bl}}{\partial t} = W_b - \frac{{\overline{w^\prime b^\prime}}_\mathrm{ent}}{\Delta B_\mathrm{bl}} 
     724\end{equation} 
     725where $h_\mathrm{bl}$ is the horizontally-varying depth of the OSBL, 
     726$\mathbf{U}_b$ and $W_b$ are the mean horizontal and vertical 
     727velocities at the base of the OSBL, ${\overline{w^\prime 
     728    b^\prime}}_\mathrm{ent}$ is the buoyancy flux due to entrainment 
     729and $\Delta B_\mathrm{bl}$ is the difference between the buoyancy 
     730averaged over the depth of the OSBL (i.e.\ including the ML and 
     731pycnocline) and the buoyancy just below the base of the OSBL. This 
     732equation for the case when the pycnocline has a finite thickness, 
     733based on the potential energy budget of the OSBL, is the leading term 
     734\citep{grant+etal18} of a generalization of that used in mixed-layer 
     735models e.g.\ \citet{kraus.turner_tellus67}, in which the thickness of the pycnocline is taken to be zero. 
     736 
     737The entrainment flux for the combination of convective and Langmuir turbulence is given by 
     738\begin{equation} \label{eq:entrain-flux} 
     739  {\overline{w^\prime b^\prime}}_\mathrm{ent} = -\alpha_{\mathrm{B}} {\overline{w^\prime b^\prime}}_0 - \alpha_{\mathrm{S}} \frac{u_*^3}{h_{\mathrm{ml}}} 
     740  + G\left(\delta/h_{\mathrm{ml}} \right)\left[\alpha_{\mathrm{S}}e^{-1.5\, \mathrm{La}_t}-\alpha_{\mathrm{L}} \frac{w_{\mathrm{*L}}^3}{h_{\mathrm{ml}}}\right] 
     741\end{equation} 
     742where the factor $G\equiv 1 - \mathrm{e}^ {-25\delta/h_{\mathrm{bl}}}(1-4\delta/h_{\mathrm{bl}})$ models the lesser efficiency of Langmuir mixing when the boundary-layer depth is much greater than the Stokes depth, and $\alpha_{\mathrm{B}}$, $\alpha_{S}$  and $\alpha_{\mathrm{L}}$ depend on the ratio of the appropriate eddy turnover time to the inertial timescale $f^{-1}$. Results from the LES suggest $\alpha_{\mathrm{B}}=0.18 F(fh_{\mathrm{bl}}/w_{*C})$, $\alpha_{S}=0.15 F(fh_{\mathrm{bl}}/u_*$  and $\alpha_{\mathrm{L}}=0.035 F(fh_{\mathrm{bl}}/u_{*L})$, where $F(x)\equiv\tanh(x^{-1})^{0.69}$. 
     743 
     744For the stable boundary layer, the equation for the depth of the OSBL is: 
     745 
     746\begin{equation}\label{eq:dhdt-stable} 
     747\max\left(\Delta B_{bl},\frac{w_{*L}^2}{h_\mathrm{bl}}\right)\frac{\partial h_\mathrm{bl}}{\partial t} = \left(0.06 + 0.52\,\frac{ h_\mathrm{bl}}{L_L}\right) \frac{w_{*L}^3}{h_\mathrm{bl}} +\left<\overline{w^\prime b^\prime}\right>_L. 
     748\end{equation} 
     749 
     750Equation. \ref{eq:dhdt-unstable} always leads to the depth of the entraining OSBL increasing (ignoring the effect of the mean vertical motion), but the change in the thickness of the stable OSBL given by Eq. \ref{eq:dhdt-stable} can be positive or negative, depending on the magnitudes of $\left<\overline{w^\prime b^\prime}\right>_L$ and $h_\mathrm{bl}/L_L$. The rate at which the depth of the OSBL can decrease is limited by choosing an effective buoyancy $w_{*L}^2/h_\mathrm{bl}$, in place of $\Delta B_{bl}$ which will be $\approx 0$ for the collapsing OSBL. 
     751 
    546752 
    547753%% ================================================================================================= 
     
    551757\begin{figure}[!t] 
    552758  \centering 
    553   \includegraphics[width=0.66\textwidth]{Fig_ZDF_TKE_time_scheme} 
     759  \includegraphics[width=0.66\textwidth]{ZDF_TKE_time_scheme} 
    554760  \caption[Subgrid kinetic energy integration in GLS and TKE schemes]{ 
    555761    Illustration of the subgrid kinetic energy integration in GLS and TKE schemes and 
     
    663869\begin{figure}[!htb] 
    664870  \centering 
    665   \includegraphics[width=0.66\textwidth]{Fig_npc} 
     871  \includegraphics[width=0.66\textwidth]{ZDF_npc} 
    666872  \caption[Unstable density profile treated by the non penetrative convective adjustment algorithm]{ 
    667873    Example of an unstable density profile treated by 
     
    8081014\begin{figure}[!t] 
    8091015  \centering 
    810   \includegraphics[width=0.66\textwidth]{Fig_zdfddm} 
     1016  \includegraphics[width=0.66\textwidth]{ZDF_ddm} 
    8111017  \caption[Diapycnal diffusivities for temperature and salt in regions of salt fingering and 
    8121018  diffusive convection]{ 
     
    12861492\begin{figure}[!t] 
    12871493  \centering 
    1288   \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_coeff} 
     1494  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_coeff} 
    12891495  \caption[Partitioning coefficient used to partition vertical velocities into parts]{ 
    12901496    The value of the partitioning coefficient (\cf) used to partition vertical velocities into 
     
    13261532\begin{figure}[!t] 
    13271533  \centering 
    1328   \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 
     1534  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_frames} 
    13291535  \caption[OVERFLOW: time-series of temperature vertical cross-sections]{ 
    13301536    A time-series of temperature vertical cross-sections for the OVERFLOW test case. 
     
    14061612\begin{figure}[!t] 
    14071613  \centering 
    1408   \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt} 
     1614  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_overflow_all_rdt} 
    14091615  \caption[OVERFLOW: sample temperature vertical cross-sections from mid- and end-run]{ 
    14101616    Sample temperature vertical cross-sections from mid- and end-run using 
     
    14191625\begin{figure}[!t] 
    14201626  \centering 
    1421   \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 
     1627  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf} 
    14221628  \caption[OVERFLOW: maximum partitioning coefficient during a series of test runs]{ 
    14231629    The maximum partitioning coefficient during a series of test runs with 
     
    14301636\begin{figure}[!t] 
    14311637  \centering 
    1432   \includegraphics[width=0.66\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 
     1638  \includegraphics[width=0.66\textwidth]{ZDF_zad_Aimp_maxCf_loc} 
    14331639  \caption[OVERFLOW: maximum partitioning coefficient for the case overlaid]{ 
    14341640    The maximum partitioning coefficient for the \forcode{nn_rdt=10.0} case overlaid with 
     
    14371643\end{figure} 
    14381644 
    1439 \onlyinsubfile{\input{../../global/epilogue}} 
     1645\subinc{\input{../../global/epilogue}} 
    14401646 
    14411647\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_cfgs.tex

    r11598 r11954  
    9393\begin{figure}[!t] 
    9494  \centering 
    95   \includegraphics[width=0.66\textwidth]{Fig_ORCA_NH_mesh} 
     95  \includegraphics[width=0.66\textwidth]{CFGS_ORCA_NH_mesh} 
    9696  \caption[ORCA mesh conception]{ 
    9797    ORCA mesh conception. 
     
    120120\begin{figure}[!tbp] 
    121121  \centering 
    122   \includegraphics[width=0.66\textwidth]{Fig_ORCA_NH_msh05_e1_e2} 
    123   \includegraphics[width=0.66\textwidth]{Fig_ORCA_aniso} 
     122  \includegraphics[width=0.66\textwidth]{CFGS_ORCA_NH_msh05_e1_e2} 
     123  \includegraphics[width=0.66\textwidth]{CFGS_ORCA_aniso} 
    124124  \caption[Horizontal scale factors and ratio of anisotropy for ORCA 0.5\deg\ mesh]{ 
    125125    \textit{Top}: Horizontal scale factors ($e_1$, $e_2$) and 
     
    264264\begin{figure}[!t] 
    265265  \centering 
    266   \includegraphics[width=0.66\textwidth]{Fig_GYRE} 
     266  \includegraphics[width=0.66\textwidth]{CFGS_GYRE} 
    267267  \caption[Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54]{ 
    268268    Snapshot of relative vorticity at the surface of the model domain in GYRE R9, R27 and R54. 
     
    292292Unlike ordinary river points the Baltic inputs also include salinity and temperature data. 
    293293 
    294 \onlyinsubfile{\input{../../global/epilogue}} 
     294\subinc{\input{../../global/epilogue}} 
    295295 
    296296\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_conservation.tex

    r11598 r11954  
    334334It has not been implemented. 
    335335 
    336 \onlyinsubfile{\input{../../global/epilogue}} 
     336\subinc{\input{../../global/epilogue}} 
    337337 
    338338\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_misc.tex

    r11598 r11954  
    9494\begin{figure}[!tbp] 
    9595  \centering 
    96   \includegraphics[width=0.66\textwidth]{Fig_Gibraltar} 
    97   \includegraphics[width=0.66\textwidth]{Fig_Gibraltar2} 
     96  \includegraphics[width=0.66\textwidth]{MISC_Gibraltar} 
     97  \includegraphics[width=0.66\textwidth]{MISC_Gibraltar2} 
    9898  \caption[Two methods to defined the Gibraltar strait]{ 
    9999    Example of the Gibraltar strait defined in a 1\deg\ $\times$ 1\deg\ mesh. 
     
    111111\begin{figure}[!tbp] 
    112112  \centering 
    113   \includegraphics[width=0.66\textwidth]{Fig_closea_mask_example} 
     113  \includegraphics[width=0.66\textwidth]{MISC_closea_mask_example} 
    114114  \caption[Mask fields for the \protect\mdl{closea} module]{ 
    115115    Example of mask fields for the \protect\mdl{closea} module. 
     
    415415increment also applies to the time.step file which is otherwise updated every timestep. 
    416416 
    417 \onlyinsubfile{\input{../../global/epilogue}} 
     417\subinc{\input{../../global/epilogue}} 
    418418 
    419419\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_model_basics.tex

    r11608 r11954  
    1313 
    1414{\footnotesize 
    15   \begin{tabularx}{\textwidth}{l||X|X} 
    16     Release & Author(s) & Modifications \\ 
     15  \begin{tabular}{l||l|l} 
     16    Release          & Author(s)                                   & Modifications      \\ 
    1717    \hline 
    18     {\em   4.0} & {\em Mike Bell                       } & {\em Update       } \\ 
    19     {\em   3.6} & {\em Gurvan Madec                    } & {\em Minor changes} \\ 
    20     {\em <=3.4} & {\em Gurvan Madec and Steven Alderson} & {\em First version} \\ 
    21   \end{tabularx} 
     18    {\em        4.0} & {\em Mike Bell                            } & {\em Review       } \\ 
     19    {\em        3.6} & {\em Tim Graham and Gurvan Madec          } & {\em Updates      } \\ 
     20    {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ 
     21  \end{tabular} 
    2222} 
    2323 
     
    130130\begin{figure} 
    131131  \centering 
    132   \includegraphics[width=0.66\textwidth]{Fig_I_ocean_bc} 
     132  \includegraphics[width=0.66\textwidth]{MB_ocean_bc} 
    133133  \caption[Ocean boundary conditions]{ 
    134134    The ocean is bounded by two surfaces, $z = - H(i,j)$ and $z = \eta(i,j,t)$, 
     
    208208\] 
    209209Two strategies can be considered for the surface pressure term: 
    210 \begin{enumerate*}[label={(\alph*)}] 
     210\begin{enumerate*}[label=(\textit{\alph*})] 
    211211\item introduce of a new variable $\eta$, the free-surface elevation, 
    212212for which a prognostic equation can be established and solved; 
     
    327327\begin{figure} 
    328328  \centering 
    329   \includegraphics[width=0.33\textwidth]{Fig_I_earth_referential} 
     329  \includegraphics[width=0.33\textwidth]{MB_earth_referential} 
    330330  \caption[Geographical and curvilinear coordinate systems]{ 
    331331    the geographical coordinate system $(\lambda,\varphi,z)$ and the curvilinear 
     
    486486\item [Flux form of the momentum equations] 
    487487  % \label{eq:MB_dyn_flux} 
    488   \begin{multline*} 
     488  \begin{alignat*}{2} 
    489489    % \label{eq:MB_dyn_flux_u} 
    490     \pd[u]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) \\ 
    491     - \frac{1}{e_3} \pd[(w \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} 
    492   \end{multline*} 
    493   \begin{multline*} 
     490    \pd[u]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, u)]{i} + \pd[(e_1 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(w \, u)]{k} \\ 
     491    &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_u^{\vect U} + F_u^{\vect U} 
     492  \end{alignat*} 
     493  \begin{alignat*}{2} 
    494494    % \label{eq:MB_dyn_flux_v} 
    495     \pd[v]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) \\ 
    496     - \frac{1}{e_3} \pd[(w \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} 
    497   \end{multline*} 
     495    \pd[v]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2} \lt( \pd[(e_2 \, u \, v)]{i} + \pd[(e_1 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(w \, v)]{k} \\ 
     496    &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) + D_v^{\vect U} + F_v^{\vect U} 
     497  \end{alignat*} 
    498498  where $\zeta$, the relative vorticity, is given by \autoref{eq:MB_curl_Uh} and 
    499499  $p_s$, the surface pressure, is given by: 
     
    579579an explicit computation of vertical advection relative to the moving s-surfaces. 
    580580 
    581 %\gmcomment{ 
    582 %A key point here is that the $s$-coordinate depends on $(i,j)$ ==> horizontal pressure gradient... 
     581\cmtgm{A key point here is that the $s$-coordinate depends on $(i,j)$ 
     582  ==> horizontal pressure gradient...} 
    583583The generalized vertical coordinates used in ocean modelling are not orthogonal, 
    584584which contrasts with many other applications in mathematical physics. 
     
    650650  \end{gather*} 
    651651\item [Flux form of the momentum equation] 
    652   \begin{multline*} 
     652  \begin{alignat*}{2} 
    653653    % \label{eq:MB_sco_u_flux} 
    654     \frac{1}{e_3} \pd[(e_3 \, u)]{t} = + \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) \\ 
    655     - \frac{1}{e_3} \pd[(\omega \, u)]{k} - \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 
    656   \end{multline*} 
    657   \begin{multline*} 
     654    \frac{1}{e_3} \pd[(e_3 \, u)]{t} = &+ \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, v - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[(e_2 \, e_3 \, u \, u)]{i} + \pd[(e_1 \, e_3 \, v \, u)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, u)]{k} \\ 
     655    &- \frac{1}{e_1} \pd[]{i} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o} \sigma_1 + D_u^{\vect U} + F_u^{\vect U} 
     656  \end{alignat*} 
     657  \begin{alignat*}{2} 
    658658  % \label{eq:MB_sco_v_flux} 
    659     \frac{1}{e_3} \pd[(e_3 \, v)]{t} = - \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) \\ 
    660     - \frac{1}{e_3} \pd[(\omega \, v)]{k} - \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
    661   \end{multline*} 
     659    \frac{1}{e_3} \pd[(e_3 \, v)]{t} = &- \lt[ f + \frac{1}{e_1 \; e_2} \lt( v \pd[e_2]{i} - u \pd[e_1]{j} \rt) \rt] \, u - \frac{1}{e_1 \; e_2 \; e_3} \lt( \pd[( e_2 \; e_3 \, u \, v)]{i} + \pd[(e_1 \; e_3 \, v \, v)]{j} \rt) - \frac{1}{e_3} \pd[(\omega \, v)]{k} \\ 
     660    &- \frac{1}{e_2} \pd[]{j} \lt( \frac{p_s + p_h}{\rho_o} \rt) - g \frac{\rho}{\rho_o}\sigma_2 + D_v^{\vect U} + F_v^{\vect U} 
     661  \end{alignat*} 
    662662  where the relative vorticity, $\zeta$, the surface pressure gradient, 
    663663  and the hydrostatic pressure have the same expressions as in $z$-coordinates although 
     
    680680and similar expressions are used for mixing and forcing terms. 
    681681 
    682 \gmcomment{ 
     682\cmtgm{ 
    683683  \colorbox{yellow}{ to be updated $= = >$} 
    684684  Add a few works on z and zps and s and underlies the differences between all of them 
     
    692692\begin{figure} 
    693693  \centering 
    694   \includegraphics[width=0.66\textwidth]{Fig_z_zstar} 
     694  \includegraphics[width=0.66\textwidth]{MB_z_zstar} 
    695695  \caption[Curvilinear z-coordinate systems (\{non-\}linear free-surface cases and re-scaled \zstar)]{ 
    696     \begin{enumerate*}[label={(\alph*)}] 
     696    \begin{enumerate*}[label=(\textit{\alph*})] 
    697697    \item $z$-coordinate in linear free-surface case; 
    698698    \item $z$-coordinate in non-linear free surface case; 
     
    10511051\begin{equation} 
    10521052  \label{eq:MB_iso_slopes} 
    1053   r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
    1054   r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
     1053  r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
     1054  r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
    10551055\end{equation} 
    10561056while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 
     
    10671067where $ \vect U^\ast = \lt( u^\ast,v^\ast,w^\ast \rt)$ is a non-divergent, 
    10681068eddy-induced transport velocity. This velocity field is defined by: 
    1069 \begin{gather*} 
     1069\[ 
    10701070  % \label{eq:MB_eiv} 
    10711071  u^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_1 \rt) \quad 
    1072   v^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \\ 
     1072  v^\ast =   \frac{1}{e_3}            \pd[]{k} \lt( A^{eiv} \;        \tilde{r}_2 \rt) \quad 
    10731073  w^\ast = - \frac{1}{e_1 e_2} \lt[   \pd[]{i} \lt( A^{eiv} \; e_2 \, \tilde{r}_1 \rt) 
    10741074                                     + \pd[]{j} \lt( A^{eiv} \; e_1 \, \tilde{r}_2 \rt) \rt] 
    1075 \end{gather*} 
     1075\] 
    10761076where $A^{eiv}$ is the eddy induced velocity coefficient 
    10771077(or equivalently the isoneutral thickness diffusivity coefficient), 
     
    11301130the $u$ and $v$-fields are considered as independent scalar fields, 
    11311131so that the diffusive operator is given by: 
    1132 \begin{gather*} 
     1132\[ 
    11331133  % \label{eq:MB_lapU_iso} 
    1134     D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \\ 
     1134    D_u^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla u \rt) \quad 
    11351135    D_v^{l \vect U} = \nabla . \lt( A^{lm} \; \Re \; \nabla v \rt) 
    1136 \end{gather*} 
     1136\] 
    11371137where $\Re$ is given by \autoref{eq:MB_iso_tensor}. 
    11381138It is the same expression as those used for diffusive operator on tracers. 
     
    11501150Nevertheless it is currently not available in the iso-neutral case. 
    11511151 
    1152 \onlyinsubfile{\input{../../global/epilogue}} 
     1152\subinc{\input{../../global/epilogue}} 
    11531153 
    11541154\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex

    r11599 r11954  
    147147\begin{figure}[!t] 
    148148  \centering 
    149   \includegraphics[width=0.66\textwidth]{Fig_DYN_dynspg_ts} 
     149  %\includegraphics[width=0.66\textwidth]{MBZ_DYN_dynspg_ts} 
    150150  \caption[Schematic of the split-explicit time stepping scheme for 
    151151  the barotropic and baroclinic modes, after \citet{Griffies2004?}]{ 
     
    311311In particular, this means that in filtered case, the matrix to be inverted has to be recomputed at each time-step. 
    312312 
    313 \onlyinsubfile{\input{../../global/epilogue}} 
     313\subinc{\input{../../global/epilogue}} 
    314314 
    315315\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_time_domain.tex

    r11599 r11954  
    1313 
    1414{\footnotesize 
    15   \begin{tabularx}{\textwidth}{l||X|X} 
    16     Release & Author(s) & Modifications \\ 
     15  \begin{tabularx}{0.5\textwidth}{l||X|X} 
     16    Release          & Author(s)                                       & 
     17    Modifications                                                      \\ 
    1718    \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 ...} 
     19    {\em        4.0} & {\em J\'{e}r\^{o}me Chanut \newline Tim Graham} & 
     20    {\em Review \newline Update                                      } \\ 
     21    {\em        3.6} & {\em Christian \'{E}th\'{e}                   } & 
     22    {\em Update                                                      } \\ 
     23    {\em $\leq$ 3.4} & {\em Gurvan Madec                             } & 
     24    {\em First version                                               } \\ 
    2225  \end{tabularx} 
    2326} 
     
    2629 
    2730% Missing things: 
    28 %  - daymod: definition of the time domain (nit000, nitend and the calendar) 
    29  
    30 \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, 
    31   would help  ==> to be added} 
    32  
    33 Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, 
     31% - daymod: definition of the time domain (nit000, nitend and the calendar) 
     32 
     33\cmtgm{STEVEN :maybe a picture of the directory structure in the introduction which 
     34could be referred to here, would help  ==> to be added} 
     35 
     36Having defined the continuous equations in \autoref{chap:MB}, 
     37we need now to choose a time discretization, 
    3438a key feature of an ocean model as it exerts a strong influence on the structure of the computer code 
    3539(\ie\ on its flowchart). 
    36 In the present chapter, we provide a general description of the \NEMO\  time stepping strategy and 
     40In the present chapter, we provide a general description of the \NEMO\ time stepping strategy and 
    3741the consequences for the order in which the equations are solved. 
    3842 
     
    4751\end{equation} 
    4852where $x$ stands for $u$, $v$, $T$ or $S$; 
    49 RHS is the Right-Hand-Side of the corresponding time evolution equation; 
     53RHS is the \textbf{R}ight-\textbf{H}and-\textbf{S}ide of the corresponding time evolution equation; 
    5054$\rdt$ is the time step; 
    5155and the superscripts indicate the time at which a quantity is evaluated. 
    52 Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated. 
     56Each term of the RHS is evaluated at a specific time stepping depending on 
     57the physics with which it is associated. 
    5358 
    5459The choice of the time stepping used for this evaluation is discussed below as well as 
    5560the implications for starting or restarting a model simulation. 
    5661Note that the time stepping calculation is generally performed in a single operation. 
    57 With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in 
    58 time for each term separately. 
     62With such a complex and nonlinear system of equations it would be dangerous to 
     63let a prognostic variable evolve in time for each term separately. 
    5964 
    6065The three level scheme requires three arrays for each prognostic variable. 
     
    6267The third array, although referred to as $x_a$ (after) in the code, 
    6368is usually not the variable at the after time step; 
    64 but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 
    65 The time stepping itself is performed once at each time step where implicit vertical diffusion is computed, \ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 
     69but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) 
     70prior to time-stepping the equation. 
     71The time stepping itself is performed once at each time step where 
     72implicit vertical diffusion is computed, 
     73\ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 
    6674 
    6775%% ================================================================================================= 
     
    6977\label{sec:TD_leap_frog} 
    7078 
    71 The time stepping used for processes other than diffusion is the well-known leapfrog scheme 
    72 \citep{mesinger.arakawa_bk76}. 
     79The time stepping used for processes other than diffusion is 
     80the well-known \textbf{L}eap\textbf{F}rog (LF) scheme \citep{mesinger.arakawa_bk76}. 
    7381This scheme is widely used for advection processes in low-viscosity fluids. 
    74 It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. 
     82It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at 
     83time step $t$, the now time step. 
    7584It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms, 
    7685but not for diffusion terms. 
    7786It is an efficient method that achieves second-order accuracy with 
    7887just one right hand side evaluation per time step. 
    79 Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by 
    80 amplifying the oscillations. 
     88Moreover, it does not artificially damp linear oscillatory motion 
     89nor does it produce instability by amplifying the oscillations. 
    8190These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme, 
    82 and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes. 
     91and the unsuitability of leapfrog differencing for the representation of diffusion and 
     92Rayleigh damping processes. 
    8393However, the scheme allows the coexistence of a numerical and a physical mode due to 
    8494its leading third order dispersive error. 
    8595In other words a divergence of odd and even time steps may occur. 
    86 To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter 
    87 (hereafter the LF-RA scheme). 
    88 This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 
     96To prevent it, the leapfrog scheme is often used in association with 
     97a \textbf{R}obert-\textbf{A}sselin time filter (hereafter the LF-RA scheme). 
     98This filter, 
     99first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 
    89100is a kind of laplacian diffusion in time that mixes odd and even time steps: 
    90101\begin{equation} 
     
    99110However, the second order truncation error is proportional to $\gamma$, which is small compared to 1. 
    100111Therefore, the LF-RA is a quasi second order accurate scheme. 
    101 The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, 
    102 because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. 
    103 When used with the 2nd order space centred discretisation of the advection terms in 
     112The LF-RA scheme is preferred to other time differencing schemes such as 
     113predictor corrector or trapezoidal schemes, because the user has an explicit and simple control of 
     114the magnitude of the time diffusion of the scheme. 
     115When used with the 2$^nd$ order space centred discretisation of the advection terms in 
    104116the momentum and tracer equations, LF-RA avoids implicit numerical diffusion: 
    105 diffusion is set explicitly by the user through the Robert-Asselin 
    106 filter parameter and the viscosity and diffusion coefficients. 
     117diffusion is set explicitly by the user through the Robert-Asselin filter parameter and 
     118the viscosity and diffusion coefficients. 
    107119 
    108120%% ================================================================================================= 
     
    110122\label{sec:TD_forward_imp} 
    111123 
    112 The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. 
     124The leapfrog differencing scheme is unsuitable for 
     125the representation of diffusion and damping processes. 
    113126For a tendency $D_x$, representing a diffusion term or a restoring term to a tracer climatology 
    114127(when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : 
     
    119132 
    120133This is diffusive in time and conditionally stable. 
    121 The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: 
     134The conditions for stability of second and fourth order horizontal diffusion schemes are 
     135\citep{griffies_bk04}: 
    122136\begin{equation} 
    123137  \label{eq:TD_euler_stability} 
     
    128142  \end{cases} 
    129143\end{equation} 
    130 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. 
     144where $e$ is the smallest grid size in the two horizontal directions and 
     145$A^h$ is the mixing coefficient. 
    131146The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient. 
    132147If it is not satisfied, even mildly, then the model soon becomes wildly unstable. 
    133 The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. 
     148The instability can be removed by either reducing the length of the time steps or 
     149reducing the mixing coefficient. 
    134150 
    135151For the vertical diffusion terms, a forward time differencing scheme can be used, 
    136 but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a 
    137 backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: 
     152but usually the numerical stability condition imposes a strong constraint on the time step. 
     153To overcome the stability constraint, a backward (or implicit) time differencing scheme is used. 
     154This scheme is unconditionally stable but diffusive and can be written as follows: 
    138155\begin{equation} 
    139156  \label{eq:TD_imp} 
     
    141158\end{equation} 
    142159 
    143 %%gm 
    144 %%gm   UPDATE the next paragraphs with time varying thickness ... 
    145 %%gm 
    146  
    147 This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: 
     160\cmtgm{UPDATE the next paragraphs with time varying thickness ...} 
     161 
     162This scheme is rather time consuming since it requires a matrix inversion. 
     163For example, the finite difference approximation of the temperature equation is: 
    148164\[ 
    149165  % \label{eq:TD_imp_zdf} 
     
    159175\end{equation} 
    160176where 
    161 \begin{align*} 
    162   c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k)     \\ 
    163   d(k) &= e_{3t}   (k)       \, / \, (2 \rdt) + c_k + c_{k + 1}    \\ 
    164   b(k) &= e_{3t}   (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) 
    165 \end{align*} 
    166  
    167 \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. 
    168 Moreover, 
    169 $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, 
     177\[ 
     178  c(k) = A_w^{vT} (k) \, / \, e_{3w} (k) \text{,} \quad 
     179  d(k) = e_{3t}   (k)       \, / \, (2 \rdt) + c_k + c_{k + 1} \quad \text{and} \quad 
     180  b(k) = e_{3t}   (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) 
     181\] 
     182 
     183\autoref{eq:TD_imp_mat} is a linear system of equations with 
     184an associated matrix which is tridiagonal. 
     185Moreover, $c(k)$ and $d(k)$ are positive and 
     186the diagonal term is greater than the sum of the two extra-diagonal terms, 
    170187therefore a special adaptation of the Gauss elimination procedure is used to find the solution 
    171188(see for example \citet{richtmyer.morton_bk67}). 
     
    175192\label{sec:TD_spg_ts} 
    176193 
    177 The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 
    178 at \textit{now} time step. This refers to as the explicit free surface case in the code (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 
    179 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation 
    180 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface 
    181 (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc 
    182 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of 
    183 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 
    184  
    185 Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), the use of a split-explicit free surface is advantageous 
    186 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication 
    187 time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. 
    188  
    189 %\gmcomment{ 
    190 \begin{figure}[!t] 
     194The leapfrog environment supports a centred in time computation of the surface pressure, 
     195\ie\ evaluated at \textit{now} time step. 
     196This refers to as the explicit free surface case in the code 
     197(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 
     198This choice however imposes a strong constraint on the time step which 
     199should be small enough to resolve the propagation of external gravity waves. 
     200As a matter of fact, one rather use in a realistic setup, 
     201a split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which 
     202barotropic and baroclinic dynamical equations are solved separately with ad-hoc time steps. 
     203The use of the time-splitting (in combination with non-linear free surface) imposes 
     204some constraints on the design of the overall flowchart, 
     205in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 
     206 
     207Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), 
     208the use of a split-explicit free surface is advantageous on massively parallel computers. 
     209Indeed, no global computations are anymore required by the elliptic solver which 
     210saves a substantial amount of communication time. 
     211Fast barotropic motions (such as tides) are also simulated with a better accuracy. 
     212 
     213%\cmtgm{ 
     214\begin{figure} 
    191215  \centering 
    192   \includegraphics[width=0.66\textwidth]{Fig_TimeStepping_flowchart_v4} 
     216  \includegraphics[width=0.66\textwidth]{TD_TimeStepping_flowchart_v4} 
    193217  \caption[Leapfrog time stepping sequence with split-explicit free surface]{ 
    194218    Sketch of the leapfrog time stepping sequence in \NEMO\ with split-explicit free surface. 
    195     The latter combined with non-linear free surface requires the dynamical tendency being 
    196     updated prior tracers tendency to ensure conservation. 
     219    The latter combined with non-linear free surface requires 
     220    the dynamical tendency being updated prior tracers tendency to ensure conservation. 
    197221    Note the use of time integrated fluxes issued from the barotropic loop in 
    198222    subsequent calculations of tracer advection and in the continuity equation. 
     
    203227 
    204228%% ================================================================================================= 
    205 \section{Modified Leapfrog -- Asselin filter scheme} 
     229\section{Modified LeapFrog -- Robert Asselin filter scheme (LF-RA)} 
    206230\label{sec:TD_mLF} 
    207231 
    208 Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to 
    209 ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 
     232Significant changes have been introduced by \cite{leclair.madec_OM09} in 
     233the LF-RA scheme in order to ensure tracer conservation and to 
     234allow the use of a much smaller value of the Asselin filter parameter. 
    210235The modifications affect both the forcing and filtering treatments in the LF-RA scheme. 
    211236 
    212 In a classical LF-RA environment, the forcing term is centred in time, 
    213 \ie\ it is time-stepped over a $2 \rdt$ period: 
     237In a classical LF-RA environment, 
     238the forcing term is centred in time, \ie\ it is time-stepped over a $2 \rdt$ period: 
    214239$x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 
    215 and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. 
     240and the time filter is given by \autoref{eq:TD_asselin} so that 
     241$Q$ is redistributed over several time step. 
    216242In the modified LF-RA environment, these two formulations have been replaced by: 
    217243\begin{gather} 
     
    222248                    - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 
    223249\end{gather} 
    224 The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) 
    225 has a significant effect: 
    226 the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. 
     250The change in the forcing formulation given by \autoref{eq:TD_forcing} 
     251(see \autoref{fig:TD_MLF_forcing}) has a significant effect: 
     252the forcing term no longer excites the divergence of odd and even time steps 
     253\citep{leclair.madec_OM09}. 
    227254% forcing seen by the model.... 
    228255This property improves the LF-RA scheme in two aspects. 
    229256First, the LF-RA can now ensure the local and global conservation of tracers. 
    230257Indeed, time filtering is no longer required on the forcing part. 
    231 The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter 
    232 (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 
     258The influence of the Asselin filter on the forcing is explicitly removed by 
     259adding a new term in the filter (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 
    233260Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 
    234261the modified formulation becomes conservative \citep{leclair.madec_OM09}. 
    235 Second, the LF-RA becomes a truly quasi -second order scheme. 
     262Second, the LF-RA becomes a truly quasi-second order scheme. 
    236263Indeed, \autoref{eq:TD_forcing} used in combination with a careful treatment of static instability 
    237264(\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) 
     
    245272even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term. 
    246273 
    247 \begin{figure}[!t] 
     274\begin{figure} 
    248275  \centering 
    249   \includegraphics[width=0.66\textwidth]{Fig_MLF_forcing} 
     276  \includegraphics[width=0.66\textwidth]{TD_MLF_forcing} 
    250277  \caption[Forcing integration methods for modified leapfrog (top and bottom)]{ 
    251278    Illustration of forcing integration methods. 
     
    271298\end{listing} 
    272299 
    273 The first time step of this three level scheme when starting from initial conditions is a forward step 
    274 (Euler time integration): 
     300The first time step of this three level scheme when starting from initial conditions is 
     301a forward step (Euler time integration): 
    275302\[ 
    276303  % \label{eq:TD_DOM_euler} 
    277304  x^1 = x^0 + \rdt \ \text{RHS}^0 
    278305\] 
    279 This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but 
     306This is done simply by keeping the leapfrog environment 
     307(\ie\ the \autoref{eq:TD} three level time stepping) but 
    280308setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 
    281309using half the value of a leapfrog time step ($2 \rdt$). 
     
    286314running the model for $2N$ time steps in one go, 
    287315or by performing two consecutive experiments of $N$ steps with a restart. 
    288 This requires saving two time levels and many auxiliary data in the restart files in machine precision. 
     316This requires saving two time levels and many auxiliary data in 
     317the restart files in machine precision. 
    289318 
    290319Note that the time step $\rdt$, is also saved in the restart file. 
    291 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step 
    292 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting 
    293 the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model 
    294 are defined through the  \nam{run}{run} namelist variables. 
    295 \gmcomment{ 
     320When restarting, if the time step has been changed, or 
     321one of the prognostic variables at \textit{before} time step is missing, 
     322an Euler time stepping scheme is imposed. 
     323A forward initial step can still be enforced by the user by 
     324setting the namelist variable \np[=0]{nn_euler}{nn\_euler}. 
     325Other options to control the time integration of the model are defined through 
     326the \nam{run}{run} namelist variables. 
     327 
     328\cmtgm{ 
    296329add here how to force the restart to contain only one time step for operational purposes 
    297330 
    298331add also the idea of writing several restart for seasonal forecast : how is it done ? 
    299332 
    300 verify that all namelist pararmeters are truly described 
     333verify that all namelist parameters are truly described 
    301334 
    302335a word on the check of restart  ..... 
    303336} 
    304337 
    305 \gmcomment{       % add a subsection here 
     338\cmtgm{       % add a subsection here 
    306339 
    307340%% ================================================================================================= 
     
    309342\label{subsec:TD_time} 
    310343 
    311 Options are defined through the  \nam{dom}{dom} namelist variables. 
     344Options are defined through the\nam{dom}{dom} namelist variables. 
    312345 \colorbox{yellow}{add here a few word on nit000 and nitend} 
    313346 
    314347 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} 
    315348 
    316 add a description of daymod, and the model calandar (leap-year and co) 
    317  
    318 }        %% end add 
    319  
    320 \gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme 
     349add a description of daymod, and the model calendar (leap-year and co) 
     350 
     351}     %% end add 
     352 
     353\cmtgm{       % add implicit in vvl case  and Crant-Nicholson scheme 
    321354 
    322355Implicit time stepping in case of variable volume thickness. 
     
    369402} 
    370403 
    371 \onlyinsubfile{\input{../../global/epilogue}} 
     404\subinc{\input{../../global/epilogue}} 
    372405 
    373406\end{document} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/SI3/namelists/namdyn_adv

    r11026 r11954  
    22&namdyn_adv     !   Ice advection 
    33!------------------------------------------------------------------------------ 
    4    ln_adv_Pra       = .false.         !  Advection scheme (Prather) 
    5    ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho) 
     4   ln_adv_Pra       = .true.         !  Advection scheme (Prather) 
     5   ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho) 
    66      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order) 
    77/ 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/TOP/subfiles/model_description.tex

    r11591 r11954  
    1414\chapter{Model Description} 
    1515\label{chap:ModDes} 
    16 \minitoc 
     16\chaptertoc 
    1717 
    1818\section{Basics} 
     
    8080\nlst{namtrc_ldf} 
    8181%------------------------------------------------------------------------------------------------------------- 
    82 In NEMO v4.0, the passive tracer diffusion has necessarily the same form as the active tracer diffusion, meaning that the numerical scheme must be the same. However the passive tracer mixing coefficient can be chosen as a multiple of the active ones by changing the value of \textit{rn\_ldf\_multi} in namelist \textit{namtrc\_ldf}. The choice of numerical scheme is then set  in the \ngn{namtra\_ldf} namelist for the dynamic described in section 5.2 of \citep{nemo_manual}. 
     82In NEMO v4.0, the passive tracer diffusion has necessarily the same form as the active tracer diffusion, meaning that the numerical scheme must be the same. However the passive tracer mixing coefficient can be chosen as a multiple of the active ones by changing the value of \textit{rn\_ldf\_multi} in namelist \textit{namtrc\_ldf}. The choice of numerical scheme is then set  in the \nam{namtra_ldf}{namtra\_ldf} namelist for the dynamic described in section 5.2 of \citep{nemo_manual}. 
    8383 
    8484 
     
    9595 
    9696The use of newtonian damping  to climatological fields or observations is also coded, sharing the same routine dans active tracers. Boolean variables are defined in the namelist\_top\_ref to select the tracers on which restoring is applied 
    97 Options are defined through the  \ngn{namtrc\_dmp} namelist variables. The restoring term is added when the namelist parameter \np{ln\_trcdmp} is set to true. The restoring coefficient is a three-dimensional array read in a file, which name is specified by the namelist variable \np{cn\_resto\_tr}. This netcdf file can be generated using the DMP\_TOOLS tool. 
     97Options are defined through the \nam{namtrc_dmp}{namtrc\_dmp} namelist variables. The restoring term is added when the namelist parameter \np{ln\_trcdmp} is set to true. The restoring coefficient is a three-dimensional array read in a file, which name is specified by the namelist variable \np{cn\_resto\_tr}. This netcdf file can be generated using the DMP\_TOOLS tool. 
    9898 
    9999 \subsection{ Tracer positivity} 
     
    104104 
    105105Sometimes, numerical scheme can generates negative values of passive tracers concentration that must be positive. For exemple,  isopycnal diffusion can created extrema. The trcrad routine artificially corrects negative concentrations with a very crude solution that either sets negative concentration to zero without adjusting the tracer budget, or by removing negative concentration and keeping mass conservation. 
    106 The treatment of negative concentrations is an option and can be selected in the namelist \ngn{namtrc\_rad} by setting the parameter \np{ln\_trcrad}  to true. 
     106The treatment of negative concentrations is an option and can be selected in the namelist \nam{namtrc_rad}{namtrc\_rad} by setting the parameter \np{ln\_trcrad}  to true. 
    107107 
    108108\section{The SMS modules} 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/document.tex

    r11591 r11954  
    1818%% End of common preamble between main and sub-files 
    1919%% Override custom cmds for full manual compilation 
    20 \newcommand{\onlyinsubfile}[1]{#1} 
    21 \newcommand{\notinsubfile}[1]{} 
     20\newcommand{\subinc}[1]{#1} 
     21\newcommand{\subexc}[1]{} 
    2222 
    2323\begin{document} 
    2424 
    25 \renewcommand{\onlyinsubfile}[1]{} 
    26 \renewcommand{\notinsubfile}[1]{#1} 
     25\renewcommand{\subinc}[1]{} 
     26\renewcommand{\subexc}[1]{#1} 
    2727 
    2828 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/frontpage.tex

    r11591 r11954  
    3333    \LARGE Version \version\ -\ \today \\ 
    3434    \medskip 
    35     \href{http://doi.org/10.5281/zenodo.\zid}{ \includegraphics{{badges/zenodo.\zid}.pdf} } 
     35    \href{http://doi.org/10.5281/zenodo.\zid}{ \includegraphics{badges/zenodo.\zid} } 
    3636  \end{center} 
    3737 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/new_cmds.tex

    r11584 r11954  
    3434 
    3535%% Gurvan's comments 
    36 \newcommand{\gmcomment}[1]{} 
     36\newcommand{\cmtgm}[1]{} 
    3737 
    3838%% Maths 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/packages.tex

    r11595 r11954  
    1818%% Issue with fontawesome pkg: path to FontAwesome.otf has to be hard-coded 
    1919\defaultfontfeatures{ 
    20     Path = /home/ntmlod/.local/texlive2019/texmf-dist/fonts/opentype/public/fontawesome/ 
     20    Path = /usr/local/texlive/2019/texmf-dist/fonts/opentype/public/fontawesome/ 
    2121} 
    2222\usepackage{academicons, fontawesome, newtxtext} 
    2323 
    2424%% Formatting 
     25\usepackage[inline]{enumitem} 
    2526\usepackage{etoc, tabularx, xcolor} 
    2627 
    2728%% Graphics 
    28 \usepackage{caption, graphicx} 
     29\usepackage{caption, graphicx, grffile} 
    2930 
    3031%% Labels 
     
    3233 
    3334%% Mathematics 
    34 \usepackage{amsmath, amssymb} 
     35\usepackage{amsmath, amssymb, mathtools} 
    3536 
    3637%% Versatility 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/styles.tex

    r11572 r11954  
    7272} 
    7373 
     74%% Temporary fix 
     75\def\set@curr@file#1{% 
     76  \begingroup 
     77    \escapechar\m@ne 
     78    \xdef\@curr@file{\expandafter\string\csname #1\endcsname}% 
     79  \endgroup 
     80} 
     81\def\quote@name#1{"\quote@@name#1\@gobble""} 
     82\def\quote@@name#1"{#1\quote@@name} 
     83\def\unquote@name#1{\quote@@name#1\@gobble"} 
     84 
    7485\makeatother 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/nam_diadct

    r11536 r11954  
    11!----------------------------------------------------------------------- 
    2 &nam_diadct     !   transports through some sections                    (default: OFF) 
     2&nam_diadct    !   transports through some sections                     (default: OFF) 
    33!----------------------------------------------------------------------- 
    44    ln_diadct  = .false.   ! Calculate transport thru sections or not 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/nam_diaharm

    r11536 r11954  
    77       nstep_han  = 15     !    Time step frequency for harmonic analysis 
    88       tname(1)   = 'M2'   !    Name of tidal constituents 
    9        tname(2)   = 'K1' 
     9       tname(2)   = 'K1'   !              --- 
    1010/ 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/nambdy

    r11536 r11954  
    2626   ln_dyn3d_dmp  =.false.     !  open boundary condition for baroclinic velocities 
    2727   rn_time_dmp   =  1.        !  Damping time scale in days 
    28    rn_time_dmp_out = 1.        !  Outflow damping time scale 
     28   rn_time_dmp_out = 1.       !  Outflow damping time scale 
    2929   nn_rimwidth   = 10         !  width of the relaxation zone 
    3030   ln_vol        = .false.    !  total volume correction (see nn_volctl parameter) 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/nambdy_dta

    r11536 r11954  
    1111   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    1212   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    13    bn_ssh      = 'amm12_bdyT_u2d'        ,         24        , 'sossheig',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    14    bn_u2d      = 'amm12_bdyU_u2d'        ,         24        , 'vobtcrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    15    bn_v2d      = 'amm12_bdyV_u2d'        ,         24        , 'vobtcrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    16    bn_u3d      = 'amm12_bdyU_u3d'        ,         24        , 'vozocrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    17    bn_v3d      = 'amm12_bdyV_u3d'        ,         24        , 'vomecrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    18    bn_tem      = 'amm12_bdyT_tra'        ,         24        , 'votemper',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    19    bn_sal      = 'amm12_bdyT_tra'        ,         24        , 'vosaline',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     13   bn_ssh      = 'amm12_bdyT_u2d'        ,         24.       , 'sossheig',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     14   bn_u2d      = 'amm12_bdyU_u2d'        ,         24.       , 'vobtcrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     15   bn_v2d      = 'amm12_bdyV_u2d'        ,         24.       , 'vobtcrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     16   bn_u3d      = 'amm12_bdyU_u3d'        ,         24.       , 'vozocrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     17   bn_v3d      = 'amm12_bdyV_u3d'        ,         24.       , 'vomecrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     18   bn_tem      = 'amm12_bdyT_tra'        ,         24.       , 'votemper',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     19   bn_sal      = 'amm12_bdyT_tra'        ,         24.       , 'vosaline',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    2020!* for si3 
    2121   bn_a_i      = 'amm12_bdyT_ice'        ,         24.       , 'siconc'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/namberg

    r11005 r11954  
    3737   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    3838   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    39    sn_icb     =  'calving'              ,         -1         ,'calvingmask',  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
     39   sn_icb     =  'calving'              ,         -1.        ,'calvingmask',  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    4040/ 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/namc1d_uvd

    r10075 r11954  
    1010   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    1111   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    12    sn_ucur     = 'ucurrent'              ,         -1        ,'u_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Ume'   , '' 
    13    sn_vcur     = 'vcurrent'              ,         -1        ,'v_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Vme'   , '' 
     12   sn_ucur     = 'ucurrent'              ,         -1.       ,'u_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Ume'   , '' 
     13   sn_vcur     = 'vcurrent'              ,         -1.       ,'v_current',   .false.   , .true. , 'monthly' ,  ''              ,  'Vme'   , '' 
    1414/ 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/namdia

    r11536 r11954  
    22&namdia         !   Diagnostics 
    33!------------------------------------------------------------------------------ 
    4    ln_icediachk     = .false.         !  check online the heat, mass & salt budgets at each time step 
    5       !                               !     rate of ice spuriously gained/lost. For ex., rn_icechk=1. <=> 1mm/year, rn_icechk=0.1 <=> 1mm/10years 
    6       rn_icechk_cel =  1.             !     check at any gridcell           => stops the code if violated (and writes a file) 
    7       rn_icechk_glo =  0.1            !     check over the entire ice cover => only prints warnings 
     4   ln_icediachk     = .false.         !  check online heat, mass & salt budgets 
     5      !                               !   rate of ice spuriously gained/lost at each time step => rn_icechk=1 <=> 1.e-6 m/hour 
     6      rn_icechk_cel =  100.           !     check at each gridcell          (1.e-4m/h)=> stops the code if violated (and writes a file) 
     7      rn_icechk_glo =  1.             !     check over the entire ice cover (1.e-6m/h)=> only prints warnings 
    88   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    99   ln_icectl        = .false.         !  ice points output for debug (T or F) 
  • NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/namdta_dyn

    r11025 r11954  
    1010   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 
    1111   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! 
    12    sn_tem      = 'dyna_grid_T'           ,       120         , 'votemper'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    13    sn_sal      = 'dyna_grid_T'           ,       120         , 'vosaline'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    14    sn_mld      = 'dyna_grid_T'           ,       120         , 'somixhgt'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    15    sn_emp      = 'dyna_grid_T'           ,       120         , 'sowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    16    sn_fmf      = 'dyna_grid_T'           ,       120         , 'iowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    17    sn_ice      = 'dyna_grid_T'           ,       120         , 'soicecov'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    18    sn_qsr      = 'dyna_grid_T'           ,       120         , 'soshfldo'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    19    sn_wnd      = 'dyna_grid_T'           ,       120         , 'sowindsp'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    20    sn_uwd      = 'dyna_grid_U'           ,       120         , 'uocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    21    sn_vwd      = 'dyna_grid_V'           ,       120         , 'vocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    22    sn_wwd      = 'dyna_grid_W'           ,       120         , 'wocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    23    sn_avt      = 'dyna_grid_W'           ,       120         , 'voddmavs'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' 
    24    sn_ubl      = 'dyna_grid_U'           ,    &nbs