# Changeset 11954 for NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc

Ignore:
Timestamp:
2019-11-22T17:15:18+01:00 (18 months ago)
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
Removed

 r11154 ************************** Building the documentation ************************** .. todo:: :file:latex    : LaTeX sources and Latexmk configuration to build reference manuals with :file:manual_build.sh :file:namelists: Namelist blocks included in the documentation :file:rst      : |RST man|_ sources and Sphinx configuration to build this guide hereby with :file:guide_build.sh .. |RST man| replace:: reStructuredText (rst) .. warning::
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/main/bibliography.bib

 r11336 } @article{         kraus.turner_tellus67, author = {Kraus, E.B. and Turner, J.}, journal = {Tellus}, pages = {98--106}, title = {A one dimensional model of the seasonal thermocline {II}. {T}he general theory and its consequences}, volume = {19}, year = {1967} } @article{        large.ea_RG97, author        = "Large, W. G. and McWilliams, J. C. and Doney, S. C.", doi           = "10.1029/94RG01872", journal       = "Reviews of Geophysics", number        = {4}, pages         = {363--403}, publisher     = {AGU}, title         = "Oceanic vertical mixing: {A} review and a model with a nonlocal boundary layer parameterization", year          = "1994" } @techreport{      large.yeager_rpt04, title         = "Diurnal to decadal global forcing for ocean and sea-ice } @article{mcwilliams.ea_JFM97, author = {McWilliams, James C. and Sullivan, Peter P. and Moeng, Chin-Hoh}, doi = {10.1017/S0022112096004375}, journal = {Journal of Fluid Mechanics}, pages = {1--30}, title = {Langmuir turbulence in the ocean}, volume = {334}, year = {1997}, } @article{         mellor.blumberg_JPO04, 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 \begin{figure}[!tb] \centering \includegraphics[width=0.66\textwidth]{Fig_zgr} \includegraphics[width=0.66\textwidth]{DOMCFG_zgr} \caption[DOMAINcfg: default vertical mesh for ORCA2]{ Default vertical mesh for ORCA2: 30 ocean levels (L30). \begin{figure}[!ht] \centering \includegraphics[width=0.66\textwidth]{Fig_sco_function} \includegraphics[width=0.66\textwidth]{DOMCFG_sco_function} \caption[DOMAINcfg: examples of the stretching function applied to a seamount]{ Examples of the stretching function applied to a seamount; \begin{figure}[!ht] \centering \includegraphics[width=0.66\textwidth]{Fig_DOM_compare_coordinates_surface} \includegraphics[width=0.66\textwidth]{DOMCFG_compare_coordinates_surface} \caption[DOMAINcfg: comparison of $s$- and $z$-coordinate]{ A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines), This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO\ web site. \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_algos.tex

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

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

 r11599 \clearpage %%%  Appendix put in gmcomment as it has not been updated for \zstar and s coordinate %%%  Appendix put in cmtgm as it has not been updated for \zstar and s coordinate %I'm writting this appendix. It will be available in a forthcoming release of the documentation %\gmcomment{ %\cmtgm{ %% ================================================================================================= %gm comment \gmcomment{ \cmtgm{ The last equality comes from the following equation, \begin{flalign*} \label{subsec:INVARIANTS_2.6} \gmcomment{ \cmtgm{ A pressure gradient has no contribution to the evolution of the vorticity as the curl of a gradient is zero. In the $z$-coordinate, this property is satisfied locally on a C-grid with 2nd order finite differences %gm comment \gmcomment{ \cmtgm{ \begin{flalign*} \sum\limits_{i,j,k} \biggl\{   p_t\;\partial_t b_t   \biggr\}                                &&&\\ %} \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/apdx_s_coord.tex

 r11599 the expression of the 3D divergence in the $s-$coordinates established above. \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}

 r11598 \begin{figure}[tb] \centering \includegraphics[width=0.66\textwidth]{Fig_GRIFF_triad_fluxes} \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_triad_fluxes} \caption[Triads arrangement and tracer gradients to give lateral and vertical tracer fluxes]{ (a) Arrangement of triads $S_i$ and tracer gradients to \begin{figure}[tb] \centering \includegraphics[width=0.66\textwidth]{Fig_GRIFF_qcells} \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_qcells} \caption[Triad notation for quarter cells]{ Triad notation for quarter cells. \begin{figure}[h] \centering \includegraphics[width=0.66\textwidth]{Fig_GRIFF_bdry_triads} \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_bdry_triads} \caption[Boundary triads]{ (a) Uppermost model layer $k=1$ with $i,1$ and $i+1,1$ tracer points (black dots), \begin{figure}[h] \centering \includegraphics[width=0.66\textwidth]{Fig_GRIFF_MLB_triads} \includegraphics[width=0.66\textwidth]{TRIADS_GRIFF_MLB_triads} \caption[Definition of mixed-layer depth and calculation of linearly tapered triads]{ Definition of mixed-layer depth and calculation of linearly tapered triads. \] \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_ASM.tex

 r11599 \end{clines} \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DIA.tex

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

 r11599 \] \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_DOM.tex

 r11598 \label{chap:DOM} % Missing things: %  - istate: description of the initial state   ==> this has to be put elsewhere.. %                  perhaps in MISC ?  By the way the initialisation of T S and dynamics %                  should be put outside of DOM routine (better with TRC staff and off-line %                  tracers) %  -geo2ocean:  how to switch from geographic to mesh coordinate %     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled %    {\em 4.0} & {\em Simon M\"{u}ller \& Andrew Coward} & %    {\em %      Compatibility changes Major simplification has moved many of the options to external domain configuration tools. %      (see \autoref{apdx:DOMCFG}) %    }                                                                                            \\ %    {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} & %    {\em First version}                                                                          \\ % Missing things % -    istate: description of the initial state   ==> this has to be put elsewhere.. %              perhaps in MISC ?  By the way the initialisation of T S and dynamics %              should be put outside of DOM routine (better with TRC staff and off-line %              tracers) % - geo2ocean: how to switch from geographic to mesh coordinate % -    domclo: closed sea and lakes.... %              management of closea sea area: specific to global cfg, both forced and coupled \thispagestyle{plain} {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} \begin{tabularx}{0.8\textwidth}{l||X|X} Release                                                                                 & Author(s)                                                                               & Modifications                                                                           \\ \hline {\em 4.0                                                                              } & {\em Simon M\"{u}ller \& Andrew Coward \newline \newline Simona Flavoni and Tim Graham                                                       } & {\em Compatibility changes: many options moved to external domain configuration tools (see \autoref{apdx:DOMCFG}). \newline Updates                                                                             } \\ {\em 3.6                                                                              } & {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec         } & {\em Updates                                                                          } \\ {\em $\leq$ 3.4                                                                       } & {\em Gurvan Madec and S\'{e}bastien Masson                                            } & {\em First version                                                                    } \end{tabularx} } \clearpage Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, we need to choose a grid for spatial discretisation and related numerical algorithms. In the present chapter, we provide a general description of the staggered grid used in \NEMO, \label{subsec:DOM_cell} \begin{figure}[!tb] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_cell} \includegraphics[width=0.33\textwidth]{DOM_cell} \caption[Arrangement of variables in the unit cell of space domain]{ Arrangement of variables in the unit cell of space domain. $t$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and planetary vorticities are defined.} \label{fig:DOM_cell} \end{figure} The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. Special attention has been given to the homogeneity of the solution in the three spatial directions. The arrangement of variables is the same in all directions. It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). This is the generalisation to three dimensions of the well-known C'' grid in Arakawa's classification \citep{mesinger.arakawa_bk76}. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). This is the generalisation to three dimensions of the well-known C'' grid in Arakawa's classification \citep{mesinger.arakawa_bk76}. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation while the scale factors are chosen equal to their local analytical value. Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation while the scale factors are chosen equal to their local analytical value. An important point here is that the partial derivative of the scale factors must be evaluated by centred finite difference approximation, not from their analytical expression. This preserves the symmetry of the discrete set of equations and therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). This preserves the symmetry of the discrete set of equations and therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). A similar, related remark can be made about the domain size: when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). \begin{table}[!tb] when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). \begin{table} \centering \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} \hline t & $i$ & $j$ & $k$ \\ \hline u & $i + 1/2$ & $j$ & $k$ \\ \hline v & $i$ & $j + 1/2$ & $k$ \\ \hline w & $i$ & $j$ & $k + 1/2$ \\ \hline f & $i + 1/2$ & $j + 1/2$ & $k$ \\ \hline uw   & $i + 1/2$ & $j$ & $k + 1/2$ \\ \hline vw   & $i$ & $j + 1/2$ & $k + 1/2$ \\ \hline fw   & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ \begin{tabular}{|l|l|l|l|} \hline t   & $i$ & $j$ & $k$ \\ \hline u   & $i + 1/2$ & $j$ & $k$ \\ \hline v   & $i$ & $j + 1/2$ & $k$ \\ \hline w   & $i$ & $j$ & $k + 1/2$ \\ \hline f   & $i + 1/2$ & $j + 1/2$ & $k$ \\ \hline uw  & $i + 1/2$ & $j$ & $k + 1/2$ \\ \hline vw  & $i$ & $j + 1/2$ & $k + 1/2$ \\ \hline fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ \hline \end{tabular} Location of grid-points as a function of integer or integer and a half value of the column, line or level. This indexing is only used for the writing of the semi -discrete equations. This indexing is only used for the writing of the semi-discrete equations. In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface. firstly, there is no ambiguity in the scale factors appearing in the discrete equations, since they are first introduced in the continuous equations; secondly, analytical transformations encourage good practice by the definition of smoothly varying grids (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. secondly, analytical transformations encourage good practice by the definition of smoothly varying grids (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. \begin{figure}[!t] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_zgr_e3} \includegraphics[width=0.5\textwidth]{DOM_zgr_e3} \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, \label{subsec:DOM_operators} Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at the midpoint between them are: Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at the midpoint between them are: \begin{alignat*}{2} % \label{eq:DOM_di_mi} Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. These operators have the following discrete forms in the curvilinear $s$-coordinates system: $\begin{gather*} % \label{eq:DOM_grad} \nabla q \equiv \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k$ \begin{multline*} + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\ % \label{eq:DOM_lap} \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] + \frac{1}{e_{3t}} \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] \end{multline*} Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and \end{gather*} Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and its divergence defined at $t$-points: \begin{multline} \begin{multline*} % \label{eq:DOM_curl} \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k \end{multline} \begin{equation} \end{multline*} $% \label{eq:DOM_div} \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] + \frac{1}{e_{3t}} \delta_k (a_3) \end{equation} The vertical average over the whole water column is denoted by an overbar and is for a masked field q (\ie\ a quantity that is equal to zero inside solid areas):$ The vertical average over the whole water column is denoted by an overbar and is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): \label{eq:DOM_bar} where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in the direction indicated by the subscript (here $k$). $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in the direction indicated by the subscript (here $k$). In continuous form, the following properties are satisfied: \end{gather} It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at vector points $(u,v,w)$. Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie \begin{alignat}{4} It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie \begin{alignat}{5} \label{eq:DOM_di_adj} &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\ \end{alignat} In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to \label{subsec:DOM_Num_Index} \begin{figure}[!tb] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_index_hor} \includegraphics[width=0.33\textwidth]{DOM_index_hor} \caption[Horizontal integer indexing]{ Horizontal integer indexing used in the \fortran\ code. The array representation used in the \fortran\ code requires an integer indexing. However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of integer values for $t$-points only while all the other points involve integer and a half values. However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of integer values for $t$-points only while all the other points involve integer and a half values. Therefore, a specific integer indexing has been defined for points other than $t$-points (\ie\ velocity and vorticity grid-points). Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. %% ================================================================================================= \label{subsec:DOM_Num_Index_vertical} In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in the \fortran\ code is the reverse of that used in the semi -discrete equations and given in \autoref{subsec:DOM_cell}. The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in the \fortran\ code is the reverse of that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index in contrast to the indexing on the horizontal plane where the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index (compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). Since the scale factors are chosen to be strictly positive, accommodate the opposing vertical index directions in implementation and documentation. \begin{figure}[!pt] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_index_vert} \includegraphics[width=0.33\textwidth]{DOM_index_vert} \caption[Vertical integer indexing]{ Vertical integer indexing used in the \fortran\ code. Two typical methods are available to specify the spatial domain configuration; they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters and fields are read from a netCDF input file, 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}. the domain-specific parameters and fields are read from a NetCDF input file, 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}. If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, and default versions that configure the spatial domain for the GYRE reference configuration are present in the \path{./src/OCE/USR} directory. and default versions that configure the spatial domain for the GYRE reference configuration are present in the \path{./src/OCE/USR} directory. In version 4.0 there are no longer any options for reading complex bathymetries and to run similar models with and without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to overly complex code. Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when the model domain itself can be altered by runtime selections. The code previously used to perform vertical discretisation has been incorporated into an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. The next subsections summarise the parameter and fields related to the configuration of the whole model domain. These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when the model domain itself can be altered by runtime selections. The code previously used to perform vertical discretisation has been incorporated into an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. The next subsections summarise the parameter and fields related to the configuration of the whole model domain. These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. The requirements are presented in three sections: the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), \label{subsec:DOM_size} The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. %% ================================================================================================= \label{sec:DOM_hgr_fields} The explicit specification of a range of mesh-related fields are required for the definition of a configuration. The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include: \begin{clines} int    jpiglo, jpjglo, jpkglo            /* global domain sizes                                          */ int    jperio                            /* lateral global domain b.c.                                   */ double glamt, glamu, glamv, glamf        /* geographic longitude (t,u,v and f points respectively)      */ double gphit, gphiu, gphiv, gphif        /* geographic latitude                                          */ double e1t, e1u, e1v, e1f                /* horizontal scale factors                                     */ double e2t, e2u, e2v, e2f                /* horizontal scale factors                                     */ int    jpiglo, jpjglo, jpkglo     /* global domain sizes                                    */ int    jperio                     /* lateral global domain b.c.                             */ double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ double gphit, gphiu, gphiv, gphif /* geographic latitude                                    */ double e1t, e1u, e1v, e1f         /* horizontal scale factors                               */ double e2t, e2u, e2v, e2f         /* horizontal scale factors                               */ \end{clines} \begin{clines} /* Optional:                                                    */ int    ORCA, ORCA_index                  /* configuration name, configuration resolution                 */ double e1e2u, e1e2v                      /* U and V surfaces (if grid size reduction in some straits)    */ double ff_f, ff_t                        /* Coriolis parameter (if not on the sphere)                    */ /* Optional:                                                 */ int    ORCA, ORCA_index /* configuration name, configuration resolution              */ double e1e2u, e1e2v     /* U and V surfaces (if grid size reduction in some straits) */ double ff_f, ff_t       /* Coriolis parameter (if not on the sphere)                 */ \end{clines} This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for illustrated examples). 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 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 not the volume of the cells. Doing otherwise can lead to numerical instability issues. In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but in cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area computation flag: the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. \smallskip Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation flag: \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation flag: \texttt{iff} to a non-zero value to suppress their re-computation. Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points, thus no specific arrays are defined at $w$ points. Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points, thus no specific arrays are defined at $w$ points. %% ================================================================================================= \subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} \label{subsec:DOM_zgr} \begin{listing} \nlst{namdom} In the vertical, the model mesh is determined by four things: \begin{enumerate} \item the bathymetry given in meters; \item the number of levels of the model (\jp{jpk}); \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and \item the masking system, \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. \item the bathymetry given in meters; \item the number of levels of the model (\jp{jpk}); \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and \item the masking system, \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. \end{enumerate} \begin{figure}[!tb] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_z_zps_s_sps} \includegraphics[width=0.5\textwidth]{DOM_z_zps_s_sps} \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ The ocean bottom as seen by the model: (a) $z$-coordinate with full step, (b) $z$-coordinate with partial step, (c) $s$-coordinate: terrain following representation, (d) hybrid $s-z$ coordinate, (e) hybrid $s-z$ coordinate with partial step, and (f) same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} \begin{enumerate*}[label=(\textit{\alph*})] \item $z$-coordinate with full step, \item $z$-coordinate with partial step, \item $s$-coordinate: terrain following representation, \item hybrid $s-z$ coordinate, \item hybrid $s-z$ coordinate with partial step, and \item same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). \end{enumerate*} Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} \label{fig:DOM_z_zps_s_sps} \end{figure} it is not intended to be an option which can be changed in the middle of an experiment. The one exception to this statement being the choice of linear or non-linear free surface. In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. This is computationally wasteful since it uses the structures for time-varying 3D metrics for fields that (in the linear free surface case) are fixed. However, the linear free-surface is rarely used and implementing it this way means a single configuration file can support both options. By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): the coordinate follow the time-variation of the free surface so that the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface However, the linear free-surface is rarely used and implementing it this way means a single configuration file can support both options. By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): the coordinate follow the time-variation of the free surface so that the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface (in other words, the top of the ocean in not a rigid lid). Note that settings: \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 appear to be namelist options but they are no longer truly namelist options for \NEMO. \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 appear to be namelist options but they are no longer truly namelist options for \NEMO. Their value is written to and read from the domain configuration file and they should be treated as fixed parameters for a particular configuration. They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and serve both to provide a record of the choices made whilst building the configuration and to trigger appropriate code blocks within \NEMO. They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and serve both to provide a record of the choices made whilst building the configuration and to trigger appropriate code blocks within \NEMO. These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. A further choice related to vertical coordinate concerns the presence (or not) of ocean cavities beneath ice shelves within the model domain. A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, otherwise the top, wet layer of the ocean will always be at the ocean surface. This option is currently only available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied at the ocean/ice shelf interface. Within the model, the arrays describing the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. Within the model, the arrays describing the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step. \end{clines} This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in the model regardless of the choice of vertical coordinate. This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in the model regardless of the choice of vertical coordinate. With constant z-levels, e3 metrics will be uniform across each horizontal level. In the partial step case each e3 at the \jp{bottom\_level} may vary from its horizontal neighbours. And, in s-coordinates, variations can occur throughout the water column. With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that variations occur throughout the water column with displacements related to the sea surface height. With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that variations occur throughout the water column with displacements related to the sea surface height. These variations are typically much smaller than those arising from bottom fitted coordinates. The values for vertical metrics supplied in the domain configuration file can be considered as those arising from a flat sea surface with zero elevation. The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: \begin{alignat*}{2} tmask(i,j,k) &= &  & \begin{cases} 0 &\text{if $k < top\_level(i,j)$} \\ 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 0 &\text{if $k > bottom\_level(i,j)$} \end{cases} \\ umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) \end{alignat*} \begin{align*} tmask(i,j,k) &= \begin{cases} 0 &\text{if $k < top\_level(i,j)$} \\ 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 0 &\text{if $k > bottom\_level(i,j)$} \end{cases} \\ umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\ vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\ fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\ wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\ \text{with~} wmask(i,j,1) &= tmask(i,j,1) \end{align*} Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. \label{subsec:DOM_closea} When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies (\eg\ Great Lakes, Caspian sea \dots) even if the model resolution does not allow their communication with the rest of the ocean. When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if the model resolution does not allow their communication with the rest of the ocean. This is unnecessary when the ocean is forced by fixed atmospheric conditions, so these seas can be removed from the ocean domain. The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it should be noted here that a successful use of these options requires appropriate mask fields to be present in the domain configuration file. The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it should be noted here that a successful use of these options requires appropriate mask fields to be present in the domain configuration file. Among the possibilities are: \begin{clines} int    closea_mask          /* non-zero values in closed sea areas for optional masking                  */ int    closea_mask_rnf      /* non-zero values in closed sea areas with runoff locations (precip only)  */ int    closea_mask_emp      /* non-zero values in closed sea areas with runoff locations (total emp)     */ int closea_mask     /* non-zero values in closed sea areas for optional masking                */ int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)   */ \end{clines} Most of the arrays relating to a particular ocean model configuration discussed in this chapter (grid-point position, scale factors) can be saved in a file if namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to \forcode{.true.}; (grid-point position, scale factors) can be saved in a file if namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to \forcode{.true.}; the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. This is only really useful if Alternatively, all the arrays relating to a particular ocean model configuration (grid-point position, scale factors, depths and masks) can be saved in a file called \texttt{mesh\_mask} if namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. (grid-point position, scale factors, depths and masks) can be saved in a file called \texttt{mesh\_mask} if namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. This file contains additional fields that can be useful for post-processing applications. \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} \label{sec:DOM_DTA_tsd} \begin{listing} \nlst{namtsd} \begin{description} \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. In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid \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. In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. The computation is done in the \mdl{dtatsd} module. \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}. \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}. The default version sets horizontally uniform T and profiles as used in the GYRE configuration (see \autoref{sec:CFGS_gyre}). \end{description} \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}

• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_OBS.tex

 r11598 \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_OBS_avg_rec} \includegraphics[width=0.66\textwidth]{OBS_avg_rec} \caption[Observational weights with a rectangular footprint]{ Weights associated with each model grid box (blue lines and numbers) \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_OBS_avg_rad} \includegraphics[width=0.66\textwidth]{OBS_avg_rad} \caption[Observational weights with a radial footprint]{ Weights associated with each model grid box (blue lines and numbers) \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_ASM_obsdist_local} \includegraphics[width=0.66\textwidth]{OBS_obsdist_local} \caption[Observations with the geographical distribution]{ Example of the distribution of observations with \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_ASM_obsdist_global} \includegraphics[width=0.66\textwidth]{OBS_obsdist_global} \caption[Observations with the round-robin distribution]{ Example of the distribution of observations with %% ================================================================================================= \section{Standalone observation operator} \section{Standalone observation operator (\texttt{SAO})} \label{sec:OBS_sao} \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_OBS_dataplot_main} \includegraphics[width=0.66\textwidth]{OBS_dataplot_main} \caption{Main window of dataplot} \label{fig:OBS_dataplotmain} \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_OBS_dataplot_prof} \includegraphics[width=0.66\textwidth]{OBS_dataplot_prof} \caption[Profile plot from dataplot]{ Profile plot from dataplot produced by right clicking on a point in the main window} \end{figure} \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_SBC.tex

 r11599 %ENDIF %\gmcomment{  word doc of runoffs: %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. %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. %The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between. %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: \cmtgm{  word doc of runoffs: 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. 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. The depth option makes it possible to have the river water affecting just the surface layer, throughout depth, or some specified point in between. To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the \mdl{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/m^2s$ 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:} %% ================================================================================================= Two different bulk formulae are available: \begin{description} \item [{\np[=1]{nn_isfblk}{nn\_isfblk}}]: The melt rate is based on a balance between the upward ocean heat flux and the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}. \item [{\np[=2]{nn_isfblk}{nn\_isfblk}}]: The melt rate and the heat flux are based on a 3 equations formulation (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation). A complete description is available in \citet{jenkins_JGR91}. \end{description} Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{losch_JGR08}. Its thickness is defined by \np{rn_hisf_tbl}{rn\_hisf\_tbl}. 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. Then, the fluxes are spread over the same thickness (ie over one or several cells). 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. This can lead to super-cool temperature in the top cell under melting condition. 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.\\ Each melt bulk formula depends on a exchange coeficient ($\Gamma^{T,S}$) between the ocean and the ice. There are 3 different ways to compute the exchange coeficient: \begin{description} \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}. \begin{gather*} \begin{description} \item [{\np[=1]{nn_isfblk}{nn\_isfblk}}]: The melt rate is based on a balance between the upward ocean heat flux and the latent heat flux at the ice shelf base. A complete description is available in \citet{hunter_rpt06}. \item [{\np[=2]{nn_isfblk}{nn\_isfblk}}]: The melt rate and the heat flux are based on a 3 equations formulation (a heat flux budget at the ice base, a salt flux budget at the ice base and a linearised freezing point temperature equation). A complete description is available in \citet{jenkins_JGR91}. \end{description} Temperature and salinity used to compute the melt are the average temperature in the top boundary layer \citet{losch_JGR08}. Its thickness is defined by \np{rn_hisf_tbl}{rn\_hisf\_tbl}. 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. Then, the fluxes are spread over the same thickness (ie over one or several cells). 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. This can lead to super-cool temperature in the top cell under melting condition. 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.\\ Each melt bulk formula depends on a exchange coeficient ($\Gamma^{T,S}$) between the ocean and the ice. There are 3 different ways to compute the exchange coeficient: \begin{description} \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}. \begin{gather*} % \label{eq:SBC_isf_gamma_iso} \gamma^{T} = rn\_gammat0 \\ \gamma^{S} = rn\_gammas0 \end{gather*} This is the recommended formulation for ISOMIP. \item [{\np[=1]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity dependent and defined as \begin{gather*} \gamma^{T} = rn\_gammat0 \times u_{*} \\ \gamma^{S} = rn\_gammas0 \times u_{*} \end{gather*} where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters). See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application. \item [{\np[=2]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity and stability dependent and defined as: $\gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}}$ where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters), $\Gamma_{Turb}$ the contribution of the ocean stability and $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion. See \citet{holland.jenkins_JPO99} for all the details on this formulation. This formulation has not been extensively tested in \NEMO\ (not recommended). \end{description} \item [{\np[=2]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting. The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np[=3]{nn_isf}{nn\_isf}). The effective melting length (\np{sn_Leff_isf}{sn\_Leff\_isf}) is read from a file. \item [{\np[=3]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. The fwf (\np{sn_rnfisf}{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}). The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. \item [{\np[=4]{nn_isf}{nn\_isf}}]: The ice shelf cavity is opened (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed). However, the fwf is not computed but specified from file \np{sn_fwfisf}{sn\_fwfisf}). The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 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})\\ \gamma^{T} = rn\_gammat0 \\ \gamma^{S} = rn\_gammas0 \end{gather*} This is the recommended formulation for ISOMIP. \item [{\np[=1]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity dependent and defined as \begin{gather*} \gamma^{T} = rn\_gammat0 \times u_{*} \\ \gamma^{S} = rn\_gammas0 \times u_{*} \end{gather*} where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters). See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application. \item [{\np[=2]{nn_gammablk}{nn\_gammablk}}]: The salt and heat exchange coefficients are velocity and stability dependent and defined as: $\gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}}$ where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_hisf_tbl}{rn\_hisf\_tbl} meters), $\Gamma_{Turb}$ the contribution of the ocean stability and $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion. See \citet{holland.jenkins_JPO99} for all the details on this formulation. This formulation has not been extensively tested in \NEMO\ (not recommended). \end{description} \item [{\np[=2]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting. The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}) as in (\np[=3]{nn_isf}{nn\_isf}). The effective melting length (\np{sn_Leff_isf}{sn\_Leff\_isf}) is read from a file. \item [{\np[=3]{nn_isf}{nn\_isf}}]: The ice shelf cavity is not represented. The fwf (\np{sn_rnfisf}{sn\_rnfisf}) is prescribed and distributed along the ice shelf edge between the depth of the average grounding line (GL) (\np{sn_depmax_isf}{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn_depmin_isf}{sn\_depmin\_isf}). The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. \item [{\np[=4]{nn_isf}{nn\_isf}}]: The ice shelf cavity is opened (\np[=.true.]{ln_isfcav}{ln\_isfcav} needed). However, the fwf is not computed but specified from file \np{sn_fwfisf}{sn\_fwfisf}). The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$. 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}) \end{description} \begin{figure}[!t] \centering \includegraphics[width=0.66\textwidth]{Fig_SBC_isf} \includegraphics[width=0.66\textwidth]{SBC_isf} \caption[Ice shelf location and fresh water flux definition]{ Illustration of the location where the fwf is injected and \begin{figure}[!t] \centering \includegraphics[width=0.66\textwidth]{Fig_SBC_diurnal} \includegraphics[width=0.66\textwidth]{SBC_diurnal} \caption[Reconstruction of the diurnal cycle variation of short wave flux]{ Example of reconstruction of the diurnal cycle variation of short wave flux from \begin{figure}[!t] \centering \includegraphics[width=0.66\textwidth]{Fig_SBC_dcy} \includegraphics[width=0.66\textwidth]{SBC_dcy} \caption[Reconstruction of the diurnal cycle variation of short wave flux on an ORCA2 grid]{ Example of reconstruction of the diurnal cycle variation of short wave flux from % in ocean-ice models. \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_STO.tex

 r11598 The first four parameters define the stochastic part of equation of state. \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}

• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/NEMO/subfiles/chap_model_basics_zstar.tex

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

 r11599 {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \begin{tabularx}{0.5\textwidth}{l||X|X} Release          & Author(s)                                       & Modifications                                                      \\ \hline {\em   4.0} & {\em ...} & {\em ...} \\ {\em   3.6} & {\em ...} & {\em ...} \\ {\em   3.4} & {\em ...} & {\em ...} \\ {\em <=3.4} & {\em ...} & {\em ...} {\em        4.0} & {\em J\'{e}r\^{o}me Chanut \newline Tim Graham} & {\em Review \newline Update                                      } \\ {\em        3.6} & {\em Christian \'{E}th\'{e}                   } & {\em Update                                                      } \\ {\em $\leq$ 3.4} & {\em Gurvan Madec                             } & {\em First version                                               } \\ \end{tabularx} } % Missing things: %  - daymod: definition of the time domain (nit000, nitend and the calendar) \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, would help  ==> to be added} Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, % - daymod: definition of the time domain (nit000, nitend and the calendar) \cmtgm{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, would help  ==> to be added} Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, a key feature of an ocean model as it exerts a strong influence on the structure of the computer code (\ie\ on its flowchart). In the present chapter, we provide a general description of the \NEMO\  time stepping strategy and In the present chapter, we provide a general description of the \NEMO\ time stepping strategy and the consequences for the order in which the equations are solved. where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the corresponding time evolution equation; RHS is the \textbf{R}ight-\textbf{H}and-\textbf{S}ide of the corresponding time evolution equation; $\rdt$ is the time step; and the superscripts indicate the time at which a quantity is evaluated. Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated. Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated. The choice of the time stepping used for this evaluation is discussed below as well as the implications for starting or restarting a model simulation. Note that the time stepping calculation is generally performed in a single operation. With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in time for each term separately. With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in time for each term separately. The three level scheme requires three arrays for each prognostic variable. The third array, although referred to as $x_a$ (after) in the code, is usually not the variable at the after time step; but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 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. but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 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. %% ================================================================================================= \label{sec:TD_leap_frog} The time stepping used for processes other than diffusion is the well-known leapfrog scheme \citep{mesinger.arakawa_bk76}. The time stepping used for processes other than diffusion is the well-known \textbf{L}eap\textbf{F}rog (LF) scheme \citep{mesinger.arakawa_bk76}. This scheme is widely used for advection processes in low-viscosity fluids. It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms, but not for diffusion terms. It is an efficient method that achieves second-order accuracy with just one right hand side evaluation per time step. Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by amplifying the oscillations. Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by amplifying the oscillations. These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme, and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes. and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes. However, the scheme allows the coexistence of a numerical and a physical mode due to its leading third order dispersive error. In other words a divergence of odd and even time steps may occur. To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter (hereafter the LF-RA scheme). This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, To prevent it, the leapfrog scheme is often used in association with a \textbf{R}obert-\textbf{A}sselin time filter (hereafter the LF-RA scheme). This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, is a kind of laplacian diffusion in time that mixes odd and even time steps: However, the second order truncation error is proportional to $\gamma$, which is small compared to 1. Therefore, the LF-RA is a quasi second order accurate scheme. The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. When used with the 2nd order space centred discretisation of the advection terms in The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. When used with the 2$^nd$ order space centred discretisation of the advection terms in the momentum and tracer equations, LF-RA avoids implicit numerical diffusion: diffusion is set explicitly by the user through the Robert-Asselin filter parameter and the viscosity and diffusion coefficients. diffusion is set explicitly by the user through the Robert-Asselin filter parameter and the viscosity and diffusion coefficients. %% ================================================================================================= \label{sec:TD_forward_imp} The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. For a tendency $D_x$, representing a diffusion term or a restoring term to a tracer climatology (when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : This is diffusive in time and conditionally stable. The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: \label{eq:TD_euler_stability} \end{cases} where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient. If it is not satisfied, even mildly, then the model soon becomes wildly unstable. The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. For the vertical diffusion terms, a forward time differencing scheme can be used, but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: \label{eq:TD_imp} %%gm %%gm   UPDATE the next paragraphs with time varying thickness ... %%gm This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: \cmtgm{UPDATE the next paragraphs with time varying thickness ...} This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: % \label{eq:TD_imp_zdf} where \begin{align*} c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k) \\ d(k) &= e_{3t} (k) \, / \, (2 \rdt) + c_k + c_{k + 1} \\ b(k) &= e_{3t} (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) \end{align*} \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. Moreover, c(k) and d(k) are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, \[ c(k) = A_w^{vT} (k) \, / \, e_{3w} (k) \text{,} \quad d(k) = e_{3t} (k) \, / \, (2 \rdt) + c_k + c_{k + 1} \quad \text{and} \quad b(k) = e_{3t} (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt) \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, therefore a special adaptation of the Gauss elimination procedure is used to find the solution (see for example \citet{richtmyer.morton_bk67}). \label{sec:TD_spg_ts} The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 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}). This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 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 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. %\gmcomment{ \begin{figure}[!t] The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 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}). This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 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 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. %\cmtgm{ \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_TimeStepping_flowchart_v4} \includegraphics[width=0.66\textwidth]{TD_TimeStepping_flowchart_v4} \caption[Leapfrog time stepping sequence with split-explicit free surface]{ Sketch of the leapfrog time stepping sequence in \NEMO\ with split-explicit free surface. The latter combined with non-linear free surface requires the dynamical tendency being updated prior tracers tendency to ensure conservation. The latter combined with non-linear free surface requires the dynamical tendency being updated prior tracers tendency to ensure conservation. Note the use of time integrated fluxes issued from the barotropic loop in subsequent calculations of tracer advection and in the continuity equation. %% ================================================================================================= \section{Modified Leapfrog -- Asselin filter scheme} \section{Modified LeapFrog -- Robert Asselin filter scheme (LF-RA)} \label{sec:TD_mLF} Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. The modifications affect both the forcing and filtering treatments in the LF-RA scheme. In a classical LF-RA environment, the forcing term is centred in time, \ie\ it is time-stepped over a $2 \rdt$ period: In a classical LF-RA environment, the forcing term is centred in time, \ie\ it is time-stepped over a $2 \rdt$ period: $x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. In the modified LF-RA environment, these two formulations have been replaced by: \begin{gather} - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) \end{gather} The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) has a significant effect: the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) has a significant effect: the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. % forcing seen by the model.... This property improves the LF-RA scheme in two aspects. First, the LF-RA can now ensure the local and global conservation of tracers. Indeed, time filtering is no longer required on the forcing part. The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, the modified formulation becomes conservative \citep{leclair.madec_OM09}. Second, the LF-RA becomes a truly quasi -second order scheme. Second, the LF-RA becomes a truly quasi-second order scheme. Indeed, \autoref{eq:TD_forcing} used in combination with a careful treatment of static instability (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term. \begin{figure}[!t] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_MLF_forcing} \includegraphics[width=0.66\textwidth]{TD_MLF_forcing} \caption[Forcing integration methods for modified leapfrog (top and bottom)]{ Illustration of forcing integration methods. \end{listing} The first time step of this three level scheme when starting from initial conditions is a forward step (Euler time integration): The first time step of this three level scheme when starting from initial conditions is a forward step (Euler time integration): $% \label{eq:TD_DOM_euler} x^1 = x^0 + \rdt \ \text{RHS}^0$ This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and using half the value of a leapfrog time step ($2 \rdt$). running the model for $2N$ time steps in one go, or by performing two consecutive experiments of $N$ steps with a restart. This requires saving two time levels and many auxiliary data in the restart files in machine precision. This requires saving two time levels and many auxiliary data in the restart files in machine precision. Note that the time step $\rdt$, is also saved in the restart file. When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model are defined through the  \nam{run}{run} namelist variables. \gmcomment{ When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model are defined through the \nam{run}{run} namelist variables. \cmtgm{ add here how to force the restart to contain only one time step for operational purposes add also the idea of writing several restart for seasonal forecast : how is it done ? verify that all namelist pararmeters are truly described verify that all namelist parameters are truly described a word on the check of restart  ..... } \gmcomment{       % add a subsection here \cmtgm{       % add a subsection here %% ================================================================================================= \label{subsec:TD_time} Options are defined through the  \nam{dom}{dom} namelist variables. Options are defined through the\nam{dom}{dom} namelist variables. \colorbox{yellow}{add here a few word on nit000 and nitend} \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} add a description of daymod, and the model calandar (leap-year and co) }        %% end add \gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme add a description of daymod, and the model calendar (leap-year and co) }     %% end add \cmtgm{       % add implicit in vvl case  and Crant-Nicholson scheme Implicit time stepping in case of variable volume thickness. } \onlyinsubfile{\input{../../global/epilogue}} \subinc{\input{../../global/epilogue}} \end{document}

• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/document.tex

• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/frontpage.tex

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

• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/packages.tex

 r11595 %% Issue with fontawesome pkg: path to FontAwesome.otf has to be hard-coded \defaultfontfeatures{ Path = /home/ntmlod/.local/texlive2019/texmf-dist/fonts/opentype/public/fontawesome/ Path = /usr/local/texlive/2019/texmf-dist/fonts/opentype/public/fontawesome/ } \usepackage{academicons, fontawesome, newtxtext} %% Formatting \usepackage[inline]{enumitem} \usepackage{etoc, tabularx, xcolor} %% Graphics \usepackage{caption, graphicx} \usepackage{caption, graphicx, grffile} %% Labels %% Mathematics \usepackage{amsmath, amssymb} \usepackage{amsmath, amssymb, mathtools} %% Versatility
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/latex/global/styles.tex

 r11572 } %% Temporary fix \def\set@curr@file#1{% \begingroup \escapechar\m@ne \xdef\@curr@file{\expandafter\string\csname #1\endcsname}% \endgroup } \def\quote@name#1{"\quote@@name#1\@gobble""} \def\quote@@name#1"{#1\quote@@name} \def\unquote@name#1{\quote@@name#1\@gobble"} \makeatother

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

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

 r11536 ln_dyn3d_dmp  =.false.     !  open boundary condition for baroclinic velocities rn_time_dmp   =  1.        !  Damping time scale in days rn_time_dmp_out = 1.        !  Outflow damping time scale rn_time_dmp_out = 1.       !  Outflow damping time scale nn_rimwidth   = 10         !  width of the relaxation zone 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 !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! bn_ssh      = 'amm12_bdyT_u2d'        ,         24        , 'sossheig',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_u2d      = 'amm12_bdyU_u2d'        ,         24        , 'vobtcrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_v2d      = 'amm12_bdyV_u2d'        ,         24        , 'vobtcrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_u3d      = 'amm12_bdyU_u3d'        ,         24        , 'vozocrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_v3d      = 'amm12_bdyV_u3d'        ,         24        , 'vomecrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_tem      = 'amm12_bdyT_tra'        ,         24        , 'votemper',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_sal      = 'amm12_bdyT_tra'        ,         24        , 'vosaline',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_ssh      = 'amm12_bdyT_u2d'        ,         24.       , 'sossheig',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_u2d      = 'amm12_bdyU_u2d'        ,         24.       , 'vobtcrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_v2d      = 'amm12_bdyV_u2d'        ,         24.       , 'vobtcrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_u3d      = 'amm12_bdyU_u3d'        ,         24.       , 'vozocrtx',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_v3d      = 'amm12_bdyV_u3d'        ,         24.       , 'vomecrty',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_tem      = 'amm12_bdyT_tra'        ,         24.       , 'votemper',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' bn_sal      = 'amm12_bdyT_tra'        ,         24.       , 'vosaline',    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' !* for si3 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 !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! sn_icb     =  'calving'              ,         -1         ,'calvingmask',  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_icb     =  'calving'              ,         -1.        ,'calvingmask',  .true.   , .true. , 'yearly'  , ''               , ''       , '' /
• ## NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/doc/namelists/namc1d_uvd

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

 r11536 &namdia         !   Diagnostics !------------------------------------------------------------------------------ ln_icediachk     = .false.         !  check online the heat, mass & salt budgets at each time step !                               !     rate of ice spuriously gained/lost. For ex., rn_icechk=1. <=> 1mm/year, rn_icechk=0.1 <=> 1mm/10years rn_icechk_cel =  1.             !     check at any gridcell           => stops the code if violated (and writes a file) rn_icechk_glo =  0.1            !     check over the entire ice cover => only prints warnings ln_icediachk     = .false.         !  check online heat, mass & salt budgets !                               !   rate of ice spuriously gained/lost at each time step => rn_icechk=1 <=> 1.e-6 m/hour rn_icechk_cel =  100.           !     check at each gridcell          (1.e-4m/h)=> stops the code if violated (and writes a file) rn_icechk_glo =  1.             !     check over the entire ice cover (1.e-6m/h)=> only prints warnings ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 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 !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   ! sn_tem      = 'dyna_grid_T'           ,       120         , 'votemper'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_sal      = 'dyna_grid_T'           ,       120         , 'vosaline'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_mld      = 'dyna_grid_T'           ,       120         , 'somixhgt'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_emp      = 'dyna_grid_T'           ,       120         , 'sowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_fmf      = 'dyna_grid_T'           ,       120         , 'iowaflup'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_ice      = 'dyna_grid_T'           ,       120         , 'soicecov'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_qsr      = 'dyna_grid_T'           ,       120         , 'soshfldo'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_wnd      = 'dyna_grid_T'           ,       120         , 'sowindsp'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_uwd      = 'dyna_grid_U'           ,       120         , 'uocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_vwd      = 'dyna_grid_V'           ,       120         , 'vocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_wwd      = 'dyna_grid_W'           ,       120         , 'wocetr_eff',  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_avt      = 'dyna_grid_W'           ,       120         , 'voddmavs'  ,  .true.   , .true. , 'yearly'  , ''               , ''       , '' sn_ubl      = 'dyna_grid_U'           ,    &nbs