# Changeset 11622

Ignore:
Timestamp:
2019-10-01T13:10:55+02:00 (13 months ago)
Message:

Review of "Time Domain" and DOM chapters

Location:
NEMO/trunk/doc/latex/NEMO/subfiles
Files:
3 edited

### Legend:

Unmodified
 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 ...} 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]{Fig_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]{Fig_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: $+ \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k$ \begin{multline*} $% \label{eq:DOM_lap} \Delta q \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} + \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$ 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 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} \label{eq:DOM_di_adj} \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]{Fig_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]{Fig_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]{Fig_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={(\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}).
 r11608 {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \begin{tabular}{l||l|l} Release          & Author(s)                                   & Modifications      \\ \hline {\em   4.0} & {\em Mike Bell                       } & {\em Update       } \\ {\em   3.6} & {\em Gurvan Madec                    } & {\em Minor changes} \\ {\em <=3.4} & {\em Gurvan Madec and Steven Alderson} & {\em First version} \\ \end{tabularx} {\em        4.0} & {\em Mike Bell                            } & {\em Review       } \\ {\em        3.6} & {\em Tim Graham and Gurvan Madec          } & {\em Updates      } \\ {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ \end{tabular} } \label{eq:MB_iso_slopes} r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$.
 r11599 {\footnotesize \begin{tabularx}{\textwidth}{l||X|X} Release & Author(s) & Modifications \\ \begin{tabular}{l||l|l} 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 ...} \end{tabularx} {\em        4.0} & {\em J\'{e}r\^{o}me Chanut and Tim Graham} & {\em Review       } \\ {\em        3.6} & {\em Christian \'{E}th\'{e}              } & {\em Update       } \\ {\em $\leq$ 3.4} & {\em Gurvan Madec                        } & {\em First version} \\ \end{tabular} } % 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) \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, 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 This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: 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} \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, \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. 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] \begin{figure} \centering \includegraphics[width=0.66\textwidth]{Fig_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} \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. 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{ 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  ..... \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 add a description of daymod, and the model calendar (leap-year and co) }     %% end add \gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme