New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 11622 for NEMO/trunk/doc – NEMO

Changeset 11622 for NEMO/trunk/doc


Ignore:
Timestamp:
2019-10-01T13:10:55+02:00 (5 years ago)
Author:
nicolasmartin
Message:

Review of "Time Domain" and DOM chapters

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

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex

    r11598 r11622  
    66\label{chap:DOM} 
    77 
    8 % Missing things: 
    9 %  - istate: description of the initial state   ==> this has to be put elsewhere.. 
    10 %                  perhaps in MISC ?  By the way the initialisation of T S and dynamics 
    11 %                  should be put outside of DOM routine (better with TRC staff and off-line 
    12 %                  tracers) 
    13 %  -geo2ocean:  how to switch from geographic to mesh coordinate 
    14 %     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 
    15  
    16 %    {\em 4.0} & {\em Simon M\"{u}ller \& Andrew Coward} & 
    17 %    {\em 
    18 %      Compatibility changes Major simplification has moved many of the options to external domain configuration tools. 
    19 %      (see \autoref{apdx:DOMCFG}) 
    20 %    }                                                                                            \\ 
    21 %    {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} & 
    22 %    {\em First version}                                                                          \\ 
     8% Missing things 
     9% -    istate: description of the initial state   ==> this has to be put elsewhere.. 
     10%              perhaps in MISC ?  By the way the initialisation of T S and dynamics 
     11%              should be put outside of DOM routine (better with TRC staff and off-line 
     12%              tracers) 
     13% - geo2ocean: how to switch from geographic to mesh coordinate 
     14% -    domclo: closed sea and lakes.... 
     15%              management of closea sea area: specific to global cfg, both forced and coupled 
    2316 
    2417\thispagestyle{plain} 
    … …  
    3023{\footnotesize 
    3124  \begin{tabularx}{\textwidth}{l||X|X} 
    32     Release & Author(s) & Modifications \\ 
    33     \hline 
    34     {\em   4.0} & {\em ...} & {\em ...} \\ 
    35     {\em   3.6} & {\em ...} & {\em ...} \\ 
    36     {\em   3.4} & {\em ...} & {\em ...} \\ 
    37     {\em <=3.4} & {\em ...} & {\em ...} 
     25    Release                                                                                 & 
     26    Author(s)                                                                               & 
     27    Modifications                                                                           \\ 
     28    \hline 
     29    {\em 4.0                                                                              } & 
     30    {\em Simon M\"{u}ller \& Andrew Coward \newline \newline 
     31      Simona Flavoni and Tim Graham                                                       } & 
     32    {\em Compatibility changes: many options moved to external domain configuration tools 
     33      (see \autoref{apdx:DOMCFG}). \newline 
     34      Updates                                                                             } \\ 
     35    {\em 3.6                                                                              } & 
     36    {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec         } & 
     37    {\em Updates                                                                          } \\ 
     38    {\em $\leq$ 3.4                                                                       } & 
     39    {\em Gurvan Madec and S\'{e}bastien Masson                                            } & 
     40    {\em First version                                                                    } 
    3841  \end{tabularx} 
    3942} 
    … …  
    4144\clearpage 
    4245 
    43 Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, 
     46Having defined the continuous equations in \autoref{chap:MB} and 
     47chosen a time discretisation \autoref{chap:TD}, 
    4448we need to choose a grid for spatial discretisation and related numerical algorithms. 
    4549In the present chapter, we provide a general description of the staggered grid used in \NEMO, 
    … …  
    5458\label{subsec:DOM_cell} 
    5559 
    56 \begin{figure}[!tb] 
     60\begin{figure} 
    5761  \centering 
    58   \includegraphics[width=0.66\textwidth]{Fig_cell} 
     62  \includegraphics[width=0.33\textwidth]{Fig_cell} 
    5963  \caption[Arrangement of variables in the unit cell of space domain]{ 
    6064    Arrangement of variables in the unit cell of space domain. 
    6165    $t$ indicates scalar points where 
    6266    temperature, salinity, density, pressure and horizontal divergence are defined. 
    63     $(u,v,w)$ indicates vector points, 
    64     and $f$ indicates vorticity points where 
     67    $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where 
    6568    both relative and planetary vorticities are defined.} 
    6669  \label{fig:DOM_cell} 
    6770\end{figure} 
    6871 
    69 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 
    70 centred second-order finite difference approximation. 
     72The numerical techniques used to solve the Primitive Equations in this model are based on 
     73the traditional, centred second-order finite difference approximation. 
    7174Special attention has been given to the homogeneity of the solution in the three spatial directions. 
    7275The arrangement of variables is the same in all directions. 
    73 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 
    74 the centre of each face of the cells (\autoref{fig:DOM_cell}). 
    75 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 
    76 \citep{mesinger.arakawa_bk76}. 
    77 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 
    78 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 
    79  
    80 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that 
    81 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
    82 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. 
    83 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 
    84 the grid-point where the scale factors are defined. 
     76It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with 
     77vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). 
     78This is the generalisation to three dimensions of the well-known ``C'' grid in 
     79Arakawa's classification \citep{mesinger.arakawa_bk76}. 
     80The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each 
     81vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying 
     82the $\zeta$ and $f$-points. 
     83 
     84The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by 
     85the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 
     86The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on 
     87\autoref{tab:DOM_cell}. 
     88In all the following, 
     89subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where 
     90the scale factors are defined. 
    8591Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 
    8692As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 
    8793$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 
    88 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 
    89 while the scale factors are chosen equal to their local analytical value. 
     94Discrete partial derivatives are formulated by 
     95the traditional, centred second order finite difference approximation while 
     96the scale factors are chosen equal to their local analytical value. 
    9097An important point here is that the partial derivative of the scale factors must be evaluated by 
    9198centred finite difference approximation, not from their analytical expression. 
    92 This preserves the symmetry of the discrete set of equations and therefore satisfies many of 
    93 the continuous properties (see \autoref{apdx:INVARIANTS}). 
     99This preserves the symmetry of the discrete set of equations and 
     100therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). 
    94101A similar, related remark can be made about the domain size: 
    95 when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors 
    96 (see \autoref{eq:DOM_bar} in the next section). 
    97  
    98 \begin{table}[!tb] 
     102when needed, an area, volume, or the total ocean depth must be evaluated as 
     103the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). 
     104 
     105\begin{table} 
    99106  \centering 
    100   \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 
    101     \hline 
    102     t & $i      $ & $j      $ & $k      $ \\ 
    103     \hline 
    104     u & $i + 1/2$ & $j      $ & $k      $ \\ 
    105     \hline 
    106     v & $i      $ & $j + 1/2$ & $k      $ \\ 
    107     \hline 
    108     w & $i      $ & $j      $ & $k + 1/2$ \\ 
    109     \hline 
    110     f & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
    111     \hline 
    112     uw   & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
    113     \hline 
    114     vw   & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
    115     \hline 
    116     fw   & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
     107  \begin{tabular}{|l|l|l|l|} 
     108    \hline 
     109    t   & $i      $ & $j      $ & $k      $ \\ 
     110    \hline 
     111    u   & $i + 1/2$ & $j      $ & $k      $ \\ 
     112    \hline 
     113    v   & $i      $ & $j + 1/2$ & $k      $ \\ 
     114    \hline 
     115    w   & $i      $ & $j      $ & $k + 1/2$ \\ 
     116    \hline 
     117    f   & $i + 1/2$ & $j + 1/2$ & $k      $ \\ 
     118    \hline 
     119    uw  & $i + 1/2$ & $j      $ & $k + 1/2$ \\ 
     120    \hline 
     121    vw  & $i      $ & $j + 1/2$ & $k + 1/2$ \\ 
     122    \hline 
     123    fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 
    117124    \hline 
    118125  \end{tabular} 
    … …  
    120127    Location of grid-points as a function of integer or 
    121128    integer and a half value of the column, line or level. 
    122     This indexing is only used for the writing of the semi -discrete equations. 
     129    This indexing is only used for the writing of the semi-discrete equations. 
    123130    In the code, the indexing uses integer values only and 
    124131    is positive downwards in the vertical with $k=1$ at the surface. 
    … …  
    137144firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 
    138145since they are first introduced in the continuous equations; 
    139 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 
    140 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 
     146secondly, analytical transformations encourage good practice by 
     147the definition of smoothly varying grids 
     148(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) 
     149\citep{treguier.dukowicz.ea_JGR96}. 
    141150An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 
    142 \begin{figure}[!t] 
     151\begin{figure} 
    143152  \centering 
    144   \includegraphics[width=0.66\textwidth]{Fig_zgr_e3} 
     153  \includegraphics[width=0.5\textwidth]{Fig_zgr_e3} 
    145154  \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ 
    146155    Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 
    … …  
    159168\label{subsec:DOM_operators} 
    160169 
    161 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at 
    162 the midpoint between them are: 
     170Given the values of a variable $q$ at adjacent points, 
     171the differencing and averaging operators at the midpoint between them are: 
    163172\begin{alignat*}{2} 
    164173  % \label{eq:DOM_di_mi} 
    … …  
    168177 
    169178Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 
    170 Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has 
    171 its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. 
     179Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, 
     180the gradient of a variable $q$ defined at a $t$-point has 
     181its three components defined at $u$-, $v$- and $w$-points while 
     182its Laplacian is defined at the $t$-point. 
    172183These operators have the following discrete forms in the curvilinear $s$-coordinates system: 
    173184\[ 
    … …  
    177188                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 
    178189\] 
    179 \begin{multline*} 
     190\[ 
    180191  % \label{eq:DOM_lap} 
    181192  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    … …  
    184195                  + \frac{1}{e_{3t}} 
    185196                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 
    186 \end{multline*} 
    187  
    188 Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 
    189 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 
     197\] 
     198 
     199Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, 
     200a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has 
     201its three curl components defined at $vw$-, $uw$, and $f$-points, and 
    190202its divergence defined at $t$-points: 
    191 \begin{multline} 
     203\begin{multline*} 
    192204% \label{eq:DOM_curl} 
    193205  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}} 
    … …  
    200212                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2) 
    201213                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 
    202 \end{multline} 
    203 \begin{equation} 
     214\end{multline*} 
     215\[ 
    204216% \label{eq:DOM_div} 
    205217  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 
    206218                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 
    207219                              + \frac{1}{e_{3t}} \delta_k (a_3) 
    208 \end{equation} 
    209  
    210 The vertical average over the whole water column is denoted by an overbar and is for 
    211 a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 
     220\] 
     221 
     222The vertical average over the whole water column is denoted by an overbar and 
     223is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 
    212224\begin{equation} 
    213225  \label{eq:DOM_bar} 
    … …  
    215227\end{equation} 
    216228where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 
    217 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over 
    218 all grid points of the same type in the direction indicated by the subscript (here $k$). 
     229$k^b$ and $k^o$ are the bottom and surface $k$-indices, 
     230and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in 
     231the direction indicated by the subscript (here $k$). 
    219232 
    220233In continuous form, the following properties are satisfied: 
    … …  
    226239\end{gather} 
    227240 
    228 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as 
    229 the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 
     241It is straightforward to demonstrate that these properties are verified locally in discrete form as 
     242soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 
    230243vector points $(u,v,w)$. 
    231244 
    232245Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 
    233 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 
    234 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 
    235 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
     246It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and 
     247$\delta_k$) are skew-symmetric linear operators, 
     248and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and 
     249$\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 
    236250\begin{alignat}{4} 
    237251  \label{eq:DOM_di_adj} 
    … …  
    241255\end{alignat} 
    242256 
    243 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 
     257In other words, 
     258the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 
    244259$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 
    245260These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to 
    … …  
    250265\label{subsec:DOM_Num_Index} 
    251266 
    252 \begin{figure}[!tb] 
     267\begin{figure} 
    253268  \centering 
    254   \includegraphics[width=0.66\textwidth]{Fig_index_hor} 
     269  \includegraphics[width=0.33\textwidth]{Fig_index_hor} 
    255270  \caption[Horizontal integer indexing]{ 
    256271    Horizontal integer indexing used in the \fortran\ code. 
    … …  
    261276 
    262277The array representation used in the \fortran\ code requires an integer indexing. 
    263 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 
    264 integer values for $t$-points only while all the other points involve integer and a half values. 
     278However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with 
     279the use of integer values for $t$-points only while 
     280all the other points involve integer and a half values. 
    265281Therefore, a specific integer indexing has been defined for points other than $t$-points 
    266282(\ie\ velocity and vorticity grid-points). 
    267 Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. 
     283Furthermore, the direction of the vertical indexing has been reversed and 
     284the surface level set at $k = 1$. 
    268285 
    269286%% ================================================================================================= 
    … …  
    281298\label{subsec:DOM_Num_Index_vertical} 
    282299 
    283 In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in 
    284 the \fortran\ code is the reverse of that used in the semi -discrete equations and 
    285 given in \autoref{subsec:DOM_cell}. 
    286 The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below 
    287 (\autoref{fig:DOM_index_vert}). 
     300In the vertical, the chosen indexing requires special attention since 
     301the direction of the $k$-axis in the \fortran\ code is the reverse of 
     302that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 
     303The sea surface corresponds to the $w$-level $k = 1$, 
     304which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). 
    288305The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 
    289306the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 
    290307Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 
    291308(\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 
    292 in contrast to the indexing on the horizontal plane where the $t$-point has the same index as 
    293 the nearest velocity points in the positive direction of the respective horizontal axis index 
     309in contrast to the indexing on the horizontal plane where 
     310the $t$-point has the same index as the nearest velocity points in 
     311the positive direction of the respective horizontal axis index 
    294312(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 
    295313Since the scale factors are chosen to be strictly positive, 
    … …  
    298316accommodate the opposing vertical index directions in implementation and documentation. 
    299317 
    300 \begin{figure}[!pt] 
     318\begin{figure} 
    301319  \centering 
    302   \includegraphics[width=0.66\textwidth]{Fig_index_vert} 
     320  \includegraphics[width=0.33\textwidth]{Fig_index_vert} 
    303321  \caption[Vertical integer indexing]{ 
    304322    Vertical integer indexing used in the \fortran\ code. 
    … …  
    314332 
    315333Two typical methods are available to specify the spatial domain configuration; 
    316 they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. 
     334they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in 
     335namelist \nam{cfg}{cfg}. 
    317336 
    318337If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, 
    319 the domain-specific parameters and fields are read from a netCDF input file, 
    320 whose name (without its .nc suffix) can be specified as the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 
     338the domain-specific parameters and fields are read from a NetCDF input file, 
     339whose name (without its .nc suffix) can be specified as 
     340the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 
    321341 
    322342If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, 
    … …  
    324344subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 
    325345These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, 
    326 and default versions that configure the spatial domain for the GYRE reference configuration are present in 
    327 the \path{./src/OCE/USR} directory. 
     346and default versions that configure the spatial domain for the GYRE reference configuration are 
     347present in the \path{./src/OCE/USR} directory. 
    328348 
    329349In version 4.0 there are no longer any options for reading complex bathymetries and 
    … …  
    332352to run similar models with and without partial bottom boxes and/or sigma-coordinates, 
    333353supporting such choices leads to overly complex code. 
    334 Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when 
    335 the model domain itself can be altered by runtime selections. 
    336 The code previously used to perform vertical discretisation has been incorporated into an external tool 
    337 (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 
    338  
    339 The next subsections summarise the parameter and fields related to the configuration of the whole model domain. 
    340 These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code 
    341 inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 
     354Worse still is the difficulty of ensuring the model configurations intended to be identical are 
     355indeed so when the model domain itself can be altered by runtime selections. 
     356The code previously used to perform vertical discretisation has been incorporated into 
     357an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 
     358 
     359The next subsections summarise the parameter and fields related to 
     360the configuration of the whole model domain. 
     361These represent the minimum information that must be provided either via 
     362the \np{cn_domcfg}{cn\_domcfg} file or 
     363set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 
    342364The requirements are presented in three sections: 
    343365the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), 
    … …  
    348370\label{subsec:DOM_size} 
    349371 
    350 The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for 
    351 the $i$, $j$ and $k$ directions, respectively. 
    352 Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when 
    353 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 
    354 see \autoref{sec:LBC_mpp}). 
     372The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and 
     373\jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. 
     374Note, that the variables \texttt{jpi} and \texttt{jpj} refer to 
     375the size of each processor subdomain when the code is run in parallel using domain decomposition 
     376(\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). 
    355377 
    356378The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, 
    … …  
    360382 
    361383The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. 
    362 See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. 
     384See \autoref{sec:LBC_jperio} for details on the available options and 
     385the corresponding values for \jp{jperio}. 
    363386 
    364387%% ================================================================================================= 
    … …  
    370393\label{sec:DOM_hgr_fields} 
    371394 
    372 The explicit specification of a range of mesh-related fields are required for the definition of a configuration. 
     395The explicit specification of a range of mesh-related fields are required for 
     396the definition of a configuration. 
    373397These include: 
    374398 
    375399\begin{clines} 
    376 int    jpiglo, jpjglo, jpkglo            /* global domain sizes                                          */ 
    377 int    jperio                            /* lateral global domain b.c.                                   */ 
    378 double glamt, glamu, glamv, glamf        /* geographic longitude (t,u,v and f points respectively)      */ 
    379 double gphit, gphiu, gphiv, gphif        /* geographic latitude                                          */ 
    380 double e1t, e1u, e1v, e1f                /* horizontal scale factors                                     */ 
    381 double e2t, e2u, e2v, e2f                /* horizontal scale factors                                     */ 
     400int    jpiglo, jpjglo, jpkglo     /* global domain sizes                                    */ 
     401int    jperio                     /* lateral global domain b.c.                             */ 
     402double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ 
     403double gphit, gphiu, gphiv, gphif /* geographic latitude                                    */ 
     404double e1t, e1u, e1v, e1f         /* horizontal scale factors                               */ 
     405double e2t, e2u, e2v, e2f         /* horizontal scale factors                               */ 
    382406\end{clines} 
    383407 
    … …  
    393417 
    394418\begin{clines} 
    395                                          /* Optional:                                                    */ 
    396 int    ORCA, ORCA_index                  /* configuration name, configuration resolution                 */ 
    397 double e1e2u, e1e2v                      /* U and V surfaces (if grid size reduction in some straits)    */ 
    398 double ff_f, ff_t                        /* Coriolis parameter (if not on the sphere)                    */ 
     419                        /* Optional:                                                 */ 
     420int    ORCA, ORCA_index /* configuration name, configuration resolution              */ 
     421double e1e2u, e1e2v     /* U and V surfaces (if grid size reduction in some straits) */ 
     422double ff_f, ff_t       /* Coriolis parameter (if not on the sphere)                 */ 
    399423\end{clines} 
    400424 
    … …  
    403427This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points 
    404428(see \autoref{sec:MISC_strait} for illustrated examples). 
    405 The key is to reduce the faces of $T$-cell (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 
     429The key is to reduce the faces of $T$-cell 
     430(\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 
    406431not the volume of the cells. 
    407432Doing otherwise can lead to numerical instability issues. 
    408433In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but 
    409434in cases where a gridsize reduction is required, 
    410 the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or 
    411 pre-computed in \mdl{usrdef\_hgr}. 
    412 If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. 
    413 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set 
    414 the surface-area computation flag: 
     435the unaltered surface areas at $u$ and $v$ grid points 
     436(\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. 
     437If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and 
     438the internal computation is suppressed. 
     439Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should 
     440set the surface-area computation flag: 
    415441\texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 
    416442 
    417443\smallskip 
    418444Similar logic applies to the other optional fields: 
    419 \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if 
    420 the mesh is not on a sphere. 
    421 If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 
    422 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set 
    423 the Coriolis computation flag: 
     445\texttt{ff\_f} and \texttt{ff\_t} which can be used to 
     446provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. 
     447If present these fields will be read and used and 
     448the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 
     449Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should 
     450set the Coriolis computation flag: 
    424451\texttt{iff} to a non-zero value to suppress their re-computation. 
    425452 
    426 Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points, 
    427 thus no specific arrays are defined at $w$ points. 
     453Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to 
     454those of $t$ points, thus no specific arrays are defined at $w$ points. 
    428455 
    429456%% ================================================================================================= 
    430457\subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} 
    431458\label{subsec:DOM_zgr} 
     459 
    432460\begin{listing} 
    433461  \nlst{namdom} 
    … …  
    438466In the vertical, the model mesh is determined by four things: 
    439467\begin{enumerate} 
    440   \item the bathymetry given in meters; 
    441   \item the number of levels of the model (\jp{jpk}); 
    442   \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 
    443   \item the masking system, \ie\ the number of wet model levels at each 
    444 $(i,j)$ location of the horizontal grid. 
     468\item the bathymetry given in meters; 
     469\item the number of levels of the model (\jp{jpk}); 
     470\item the analytical transformation $z(i,j,k)$ and the vertical scale factors 
     471  (derivatives of the transformation); and 
     472\item the masking system, 
     473  \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. 
    445474\end{enumerate} 
    446475 
    447 \begin{figure}[!tb] 
     476\begin{figure} 
    448477  \centering 
    449   \includegraphics[width=0.66\textwidth]{Fig_z_zps_s_sps} 
     478  \includegraphics[width=0.5\textwidth]{Fig_z_zps_s_sps} 
    450479  \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ 
    451480    The ocean bottom as seen by the model: 
    452     (a) $z$-coordinate with full step, 
    453     (b) $z$-coordinate with partial step, 
    454     (c) $s$-coordinate: terrain following representation, 
    455     (d) hybrid $s-z$ coordinate, 
    456     (e) hybrid $s-z$ coordinate with partial step, and 
    457     (f) same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 
    458     Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 
     481    \begin{enumerate*}[label={(\alph*)}] 
     482    \item $z$-coordinate with full step, 
     483    \item $z$-coordinate with partial step, 
     484    \item $s$-coordinate: terrain following representation, 
     485    \item hybrid $s-z$ coordinate, 
     486    \item hybrid $s-z$ coordinate with partial step, and 
     487    \item same as (e) but in the non-linear free surface 
     488      (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 
     489  \end{enumerate*} 
     490  Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 
    459491  \label{fig:DOM_z_zps_s_sps} 
    460492\end{figure} 
    … …  
    463495it is not intended to be an option which can be changed in the middle of an experiment. 
    464496The one exception to this statement being the choice of linear or non-linear free surface. 
    465 In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. 
     497In v4.0 the linear free surface option is implemented as 
     498a special case of the non-linear free surface. 
    466499This is computationally wasteful since it uses the structures for time-varying 3D metrics 
    467500for fields that (in the linear free surface case) are fixed. 
    468 However, the linear free-surface is rarely used and implementing it this way means 
    469 a single configuration file can support both options. 
    470  
    471 By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 
    472 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 
    473 $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 
    474 When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 
    475 the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface 
     501However, the linear free-surface is rarely used and 
     502implementing it this way means a single configuration file can support both options. 
     503 
     504By default a non-linear free surface is used 
     505(\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 
     506the coordinate follow the time-variation of the free surface so that 
     507the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 
     508When a linear free surface is assumed 
     509(\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 
     510the vertical coordinates are fixed in time, but 
     511the seawater can move up and down across the $z_0$ surface 
    476512(in other words, the top of the ocean in not a rigid lid). 
    477513 
    478514Note that settings: 
    479 \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} mentioned in the following sections 
    480 appear to be namelist options but they are no longer truly namelist options for \NEMO. 
     515\np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} 
     516mentioned in the following sections appear to be namelist options but 
     517they are no longer truly namelist options for \NEMO. 
    481518Their value is written to and read from the domain configuration file and 
    482519they should be treated as fixed parameters for a particular configuration. 
    483 They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and 
    484 serve both to provide a record of the choices made whilst building the configuration and 
    485 to trigger appropriate code blocks within \NEMO. 
     520They are namelist options for the \texttt{DOMAINcfg} tool that can be used to 
     521build the configuration file and serve both to provide a record of the choices made whilst 
     522building the configuration and to trigger appropriate code blocks within \NEMO. 
    486523These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. 
    487524 
    … …  
    501538A further choice related to vertical coordinate concerns 
    502539the presence (or not) of ocean cavities beneath ice shelves within the model domain. 
    503 A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, 
     540A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that 
     541the domain contains ocean cavities, 
    504542otherwise the top, wet layer of the ocean will always be at the ocean surface. 
    505543This option is currently only available for $z$- or $zps$-coordinates. 
    506544In the latter case, partial steps are also applied at the ocean/ice shelf interface. 
    507545 
    508 Within the model, the arrays describing the grid point depths and vertical scale factors are three set of 
    509 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 
     546Within the model, 
     547the arrays describing the grid point depths and vertical scale factors are 
     548three set of three dimensional arrays $(i,j,k)$ defined at 
     549\textit{before}, \textit{now} and \textit{after} time step. 
    510550The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 
    511551They are updated at each model time step. 
    … …  
    534574\end{clines} 
    535575 
    536 This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in 
    537 the model regardless of the choice of vertical coordinate. 
     576This set of vertical metrics is sufficient to describe the initial depth and thickness of 
     577every gridcell in the model regardless of the choice of vertical coordinate. 
    538578With constant z-levels, e3 metrics will be uniform across each horizontal level. 
    539579In the partial step case each e3 at the \jp{bottom\_level} 
    … …  
    541581may vary from its horizontal neighbours. 
    542582And, in s-coordinates, variations can occur throughout the water column. 
    543 With the non-linear free-surface, all the coordinates behave more like the s-coordinate in 
    544 that variations occur throughout the water column with displacements related to the sea surface height. 
     583With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that 
     584variations occur throughout the water column with displacements related to the sea surface height. 
    545585These variations are typically much smaller than those arising from bottom fitted coordinates. 
    546586The values for vertical metrics supplied in the domain configuration file can be considered as 
    547587those arising from a flat sea surface with zero elevation. 
    548588 
    549 The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. 
     589The \jp{bottom\_level} and \jp{top\_level} 2D arrays define 
     590the \jp{bottom\_level} and top wet levels in each grid column. 
    550591Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). 
    551592With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. 
    … …  
    556597 
    557598From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: 
    558 \begin{alignat*}{2} 
    559   tmask(i,j,k) &= &  & 
    560     \begin{cases} 
    561                   0 &\text{if $                  k  <    top\_level(i,j)$} \\ 
    562                   1 &\text{if $bottom\_level(i,j) \leq k \leq   top\_level(i,j)$} \\ 
    563                   0 &\text{if $                  k  >     bottom\_level(i,j)$} 
    564     \end{cases} 
    565   \\ 
    566   umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    567   vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
    568   fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    569                &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
    570   wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
    571   \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 
    572 \end{alignat*} 
     599\begin{align*} 
     600  tmask(i,j,k) &= 
     601  \begin{cases} 
     602    0 &\text{if $                             k <    top\_level(i,j)$} \\ 
     603    1 &\text{if $     bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 
     604    0 &\text{if $k >  bottom\_level(i,j)                            $} 
     605  \end{cases} \\ 
     606  umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     607  vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\ 
     608  fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\ 
     609  wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\ 
     610  \text{with~} wmask(i,j,1) &= tmask(i,j,1) 
     611\end{align*} 
    573612 
    574613Note that, without ice shelves cavities, 
    575 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 
    576 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 
     614masks at $t-$ and $w-$points are identical with the numerical indexing used 
     615(\autoref{subsec:DOM_Num_Index}). 
     616Nevertheless, 
     617$wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 
    577618exactly in the same way as for the bottom boundary. 
    578619 
    … …  
    588629\label{subsec:DOM_closea} 
    589630 
    590 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 
    591 (\eg\ Great Lakes, Caspian sea \dots) even if the model resolution does not allow their communication with 
    592 the rest of the ocean. 
     631When a global ocean is coupled to an atmospheric model it is better to 
     632represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if 
     633the model resolution does not allow their communication with the rest of the ocean. 
    593634This is unnecessary when the ocean is forced by fixed atmospheric conditions, 
    594635so these seas can be removed from the ocean domain. 
    595 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and 
    596 to optionally decide on the fate of any freshwater imbalance over the area. 
    597 The options are explained in \autoref{sec:MISC_closea} but it should be noted here that 
    598 a successful use of these options requires appropriate mask fields to be present in the domain configuration file. 
     636The user has the option to 
     637set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to 
     638optionally decide on the fate of any freshwater imbalance over the area. 
     639The options are explained in \autoref{sec:MISC_closea} but 
     640it should be noted here that a successful use of these options requires 
     641appropriate mask fields to be present in the domain configuration file. 
    599642Among the possibilities are: 
    600643 
    601644\begin{clines} 
    602 int    closea_mask          /* non-zero values in closed sea areas for optional masking                  */ 
    603 int    closea_mask_rnf      /* non-zero values in closed sea areas with runoff locations (precip only)  */ 
    604 int    closea_mask_emp      /* non-zero values in closed sea areas with runoff locations (total emp)     */ 
     645int closea_mask     /* non-zero values in closed sea areas for optional masking                */ 
     646int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ 
     647int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)   */ 
    605648\end{clines} 
    606649 
    … …  
    610653 
    611654Most of the arrays relating to a particular ocean model configuration discussed in this chapter 
    612 (grid-point position, scale factors) 
    613 can be saved in a file if 
    614 namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to \forcode{.true.}; 
     655(grid-point position, scale factors) can be saved in a file if 
     656namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to 
     657\forcode{.true.}; 
    615658the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. 
    616659This is only really useful if 
    … …  
    619662 
    620663Alternatively, all the arrays relating to a particular ocean model configuration 
    621 (grid-point position, scale factors, depths and masks) 
    622 can be saved in a file called \texttt{mesh\_mask} if 
    623 namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. 
     664(grid-point position, scale factors, depths and masks) can be saved in 
     665a file called \texttt{mesh\_mask} if 
     666namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to 
     667\forcode{.true.}. 
    624668This file contains additional fields that can be useful for post-processing applications. 
    625669 
    … …  
    627671\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 
    628672\label{sec:DOM_DTA_tsd} 
     673 
    629674\begin{listing} 
    630675  \nlst{namtsd} 
    … …  
    638683 
    639684\begin{description} 
    640 \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on the model grid itself or on their native input data grids. 
    641   In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 
     685\item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on 
     686  the model grid itself or on their native input data grids. 
     687  In the latter case, 
     688  the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 
    642689  (see \autoref{subsec:SBC_iof}). 
    643   The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 
     690  The information relating to the input files are specified in 
     691  the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 
    644692  The computation is done in the \mdl{dtatsd} module. 
    645 \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 
     693\item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via 
     694  a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 
    646695  The default version sets horizontally uniform T and profiles as used in the GYRE configuration 
    647696  (see \autoref{sec:CFGS_gyre}). 

  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex

    r11608 r11622  
    1313 
    1414{\footnotesize 
    15   \begin{tabularx}{\textwidth}{l||X|X} 
    16     Release & Author(s) & Modifications \\ 
     15  \begin{tabular}{l||l|l} 
     16    Release          & Author(s)                                   & Modifications      \\ 
    1717    \hline 
    18     {\em   4.0} & {\em Mike Bell                       } & {\em Update       } \\ 
    19     {\em   3.6} & {\em Gurvan Madec                    } & {\em Minor changes} \\ 
    20     {\em <=3.4} & {\em Gurvan Madec and Steven Alderson} & {\em First version} \\ 
    21   \end{tabularx} 
     18    {\em        4.0} & {\em Mike Bell                            } & {\em Review       } \\ 
     19    {\em        3.6} & {\em Tim Graham and Gurvan Madec          } & {\em Updates      } \\ 
     20    {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ 
     21  \end{tabular} 
    2222} 
    2323 
    … …  
    10511051\begin{equation} 
    10521052  \label{eq:MB_iso_slopes} 
    1053   r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
    1054   r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
     1053  r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 
     1054  r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 
    10551055\end{equation} 
    10561056while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. 

  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex

    r11599 r11622  
    1313 
    1414{\footnotesize 
    15   \begin{tabularx}{\textwidth}{l||X|X} 
    16     Release & Author(s) & Modifications \\ 
     15  \begin{tabular}{l||l|l} 
     16    Release          & Author(s)                                  & Modifications      \\ 
    1717    \hline 
    18     {\em   4.0} & {\em ...} & {\em ...} \\ 
    19     {\em   3.6} & {\em ...} & {\em ...} \\ 
    20     {\em   3.4} & {\em ...} & {\em ...} \\ 
    21     {\em <=3.4} & {\em ...} & {\em ...} 
    22   \end{tabularx} 
     18    {\em        4.0} & {\em J\'{e}r\^{o}me Chanut and Tim Graham} & {\em Review       } \\ 
     19    {\em        3.6} & {\em Christian \'{E}th\'{e}              } & {\em Update       } \\ 
     20    {\em $\leq$ 3.4} & {\em Gurvan Madec                        } & {\em First version} \\ 
     21  \end{tabular} 
    2322} 
    2423 
    … …  
    2625 
    2726% Missing things: 
    28 %  - daymod: definition of the time domain (nit000, nitend and the calendar) 
    29  
    30 \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, 
    31   would help  ==> to be added} 
    32  
    33 Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, 
     27% - daymod: definition of the time domain (nit000, nitend and the calendar) 
     28 
     29\gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which 
     30could be referred to here, would help  ==> to be added} 
     31 
     32Having defined the continuous equations in \autoref{chap:MB}, 
     33we need now to choose a time discretization, 
    3434a key feature of an ocean model as it exerts a strong influence on the structure of the computer code 
    3535(\ie\ on its flowchart). 
    36 In the present chapter, we provide a general description of the \NEMO\  time stepping strategy and 
     36In the present chapter, we provide a general description of the \NEMO\ time stepping strategy and 
    3737the consequences for the order in which the equations are solved. 
    3838 
    … …  
    4747\end{equation} 
    4848where $x$ stands for $u$, $v$, $T$ or $S$; 
    49 RHS is the Right-Hand-Side of the corresponding time evolution equation; 
     49RHS is the \textbf{R}ight-\textbf{H}and-\textbf{S}ide of the corresponding time evolution equation; 
    5050$\rdt$ is the time step; 
    5151and the superscripts indicate the time at which a quantity is evaluated. 
    52 Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated. 
     52Each term of the RHS is evaluated at a specific time stepping depending on 
     53the physics with which it is associated. 
    5354 
    5455The choice of the time stepping used for this evaluation is discussed below as well as 
    5556the implications for starting or restarting a model simulation. 
    5657Note that the time stepping calculation is generally performed in a single operation. 
    57 With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in 
    58 time for each term separately. 
     58With such a complex and nonlinear system of equations it would be dangerous to 
     59let a prognostic variable evolve in time for each term separately. 
    5960 
    6061The three level scheme requires three arrays for each prognostic variable. 
    … …  
    6263The third array, although referred to as $x_a$ (after) in the code, 
    6364is usually not the variable at the after time step; 
    64 but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 
    65 The time stepping itself is performed once at each time step where implicit vertical diffusion is computed, \ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 
     65but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) 
     66prior to time-stepping the equation. 
     67The time stepping itself is performed once at each time step where 
     68implicit vertical diffusion is computed, 
     69\ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 
    6670 
    6771%% ================================================================================================= 
    … …  
    6973\label{sec:TD_leap_frog} 
    7074 
    71 The time stepping used for processes other than diffusion is the well-known leapfrog scheme 
    72 \citep{mesinger.arakawa_bk76}. 
     75The time stepping used for processes other than diffusion is 
     76the well-known \textbf{L}eap\textbf{F}rog (LF) scheme \citep{mesinger.arakawa_bk76}. 
    7377This scheme is widely used for advection processes in low-viscosity fluids. 
    74 It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. 
     78It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at 
     79time step $t$, the now time step. 
    7580It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms, 
    7681but not for diffusion terms. 
    7782It is an efficient method that achieves second-order accuracy with 
    7883just one right hand side evaluation per time step. 
    79 Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by 
    80 amplifying the oscillations. 
     84Moreover, it does not artificially damp linear oscillatory motion 
     85nor does it produce instability by amplifying the oscillations. 
    8186These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme, 
    82 and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes. 
     87and the unsuitability of leapfrog differencing for the representation of diffusion and 
     88Rayleigh damping processes. 
    8389However, the scheme allows the coexistence of a numerical and a physical mode due to 
    8490its leading third order dispersive error. 
    8591In other words a divergence of odd and even time steps may occur. 
    86 To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter 
    87 (hereafter the LF-RA scheme). 
    88 This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 
     92To prevent it, the leapfrog scheme is often used in association with 
     93a \textbf{R}obert-\textbf{A}sselin time filter (hereafter the LF-RA scheme). 
     94This filter, 
     95first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 
    8996is a kind of laplacian diffusion in time that mixes odd and even time steps: 
    9097\begin{equation} 
    … …  
    99106However, the second order truncation error is proportional to $\gamma$, which is small compared to 1. 
    100107Therefore, the LF-RA is a quasi second order accurate scheme. 
    101 The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, 
    102 because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. 
    103 When used with the 2nd order space centred discretisation of the advection terms in 
     108The LF-RA scheme is preferred to other time differencing schemes such as 
     109predictor corrector or trapezoidal schemes, because the user has an explicit and simple control of 
     110the magnitude of the time diffusion of the scheme. 
     111When used with the 2$^nd$ order space centred discretisation of the advection terms in 
    104112the momentum and tracer equations, LF-RA avoids implicit numerical diffusion: 
    105 diffusion is set explicitly by the user through the Robert-Asselin 
    106 filter parameter and the viscosity and diffusion coefficients. 
     113diffusion is set explicitly by the user through the Robert-Asselin filter parameter and 
     114the viscosity and diffusion coefficients. 
    107115 
    108116%% ================================================================================================= 
    … …  
    110118\label{sec:TD_forward_imp} 
    111119 
    112 The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. 
     120The leapfrog differencing scheme is unsuitable for 
     121the representation of diffusion and damping processes. 
    113122For a tendency $D_x$, representing a diffusion term or a restoring term to a tracer climatology 
    114123(when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : 
    … …  
    119128 
    120129This is diffusive in time and conditionally stable. 
    121 The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: 
     130The conditions for stability of second and fourth order horizontal diffusion schemes are 
     131\citep{griffies_bk04}: 
    122132\begin{equation} 
    123133  \label{eq:TD_euler_stability} 
    … …  
    128138  \end{cases} 
    129139\end{equation} 
    130 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. 
     140where $e$ is the smallest grid size in the two horizontal directions and 
     141$A^h$ is the mixing coefficient. 
    131142The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient. 
    132143If it is not satisfied, even mildly, then the model soon becomes wildly unstable. 
    133 The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. 
     144The instability can be removed by either reducing the length of the time steps or 
     145reducing the mixing coefficient. 
    134146 
    135147For the vertical diffusion terms, a forward time differencing scheme can be used, 
    136 but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a 
    137 backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: 
     148but usually the numerical stability condition imposes a strong constraint on the time step. 
     149To overcome the stability constraint, a backward (or implicit) time differencing scheme is used. 
     150This scheme is unconditionally stable but diffusive and can be written as follows: 
    138151\begin{equation} 
    139152  \label{eq:TD_imp} 
    … …  
    145158%%gm 
    146159 
    147 This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: 
     160This scheme is rather time consuming since it requires a matrix inversion. 
     161For example, the finite difference approximation of the temperature equation is: 
    148162\[ 
    149163  % \label{eq:TD_imp_zdf} 
    … …  
    165179\end{align*} 
    166180 
    167 \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. 
    168 Moreover, 
    169 $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, 
     181\autoref{eq:TD_imp_mat} is a linear system of equations with 
     182an associated matrix which is tridiagonal. 
     183Moreover, $c(k)$ and $d(k)$ are positive and 
     184the diagonal term is greater than the sum of the two extra-diagonal terms, 
    170185therefore a special adaptation of the Gauss elimination procedure is used to find the solution 
    171186(see for example \citet{richtmyer.morton_bk67}). 
    … …  
    175190\label{sec:TD_spg_ts} 
    176191 
    177 The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 
    178 at \textit{now} time step. This refers to as the explicit free surface case in the code (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 
    179 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation 
    180 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface 
    181 (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc 
    182 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of 
    183 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 
    184  
    185 Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), the use of a split-explicit free surface is advantageous 
    186 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication 
    187 time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. 
     192The leapfrog environment supports a centred in time computation of the surface pressure, 
     193\ie\ evaluated at \textit{now} time step. 
     194This refers to as the explicit free surface case in the code 
     195(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 
     196This choice however imposes a strong constraint on the time step which 
     197should be small enough to resolve the propagation of external gravity waves. 
     198As a matter of fact, one rather use in a realistic setup, 
     199a split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which 
     200barotropic and baroclinic dynamical equations are solved separately with ad-hoc time steps. 
     201The use of the time-splitting (in combination with non-linear free surface) imposes 
     202some constraints on the design of the overall flowchart, 
     203in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 
     204 
     205Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), 
     206the use of a split-explicit free surface is advantageous on massively parallel computers. 
     207Indeed, no global computations are anymore required by the elliptic solver which 
     208saves a substantial amount of communication time. 
     209Fast barotropic motions (such as tides) are also simulated with a better accuracy. 
    188210 
    189211%\gmcomment{ 
    190 \begin{figure}[!t] 
     212\begin{figure} 
    191213  \centering 
    192214  \includegraphics[width=0.66\textwidth]{Fig_TimeStepping_flowchart_v4} 
    193215  \caption[Leapfrog time stepping sequence with split-explicit free surface]{ 
    194216    Sketch of the leapfrog time stepping sequence in \NEMO\ with split-explicit free surface. 
    195     The latter combined with non-linear free surface requires the dynamical tendency being 
    196     updated prior tracers tendency to ensure conservation. 
     217    The latter combined with non-linear free surface requires 
     218    the dynamical tendency being updated prior tracers tendency to ensure conservation. 
    197219    Note the use of time integrated fluxes issued from the barotropic loop in 
    198220    subsequent calculations of tracer advection and in the continuity equation. 
    … …  
    203225 
    204226%% ================================================================================================= 
    205 \section{Modified Leapfrog -- Asselin filter scheme} 
     227\section{Modified LeapFrog -- Robert Asselin filter scheme (LF-RA)} 
    206228\label{sec:TD_mLF} 
    207229 
    208 Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to 
    209 ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 
     230Significant changes have been introduced by \cite{leclair.madec_OM09} in 
     231the LF-RA scheme in order to ensure tracer conservation and to 
     232allow the use of a much smaller value of the Asselin filter parameter. 
    210233The modifications affect both the forcing and filtering treatments in the LF-RA scheme. 
    211234 
    212 In a classical LF-RA environment, the forcing term is centred in time, 
    213 \ie\ it is time-stepped over a $2 \rdt$ period: 
     235In a classical LF-RA environment, 
     236the forcing term is centred in time, \ie\ it is time-stepped over a $2 \rdt$ period: 
    214237$x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 
    215 and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. 
     238and the time filter is given by \autoref{eq:TD_asselin} so that 
     239$Q$ is redistributed over several time step. 
    216240In the modified LF-RA environment, these two formulations have been replaced by: 
    217241\begin{gather} 
    … …  
    222246                    - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 
    223247\end{gather} 
    224 The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) 
    225 has a significant effect: 
    226 the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. 
     248The change in the forcing formulation given by \autoref{eq:TD_forcing} 
     249(see \autoref{fig:TD_MLF_forcing}) has a significant effect: 
     250the forcing term no longer excites the divergence of odd and even time steps 
     251\citep{leclair.madec_OM09}. 
    227252% forcing seen by the model.... 
    228253This property improves the LF-RA scheme in two aspects. 
    229254First, the LF-RA can now ensure the local and global conservation of tracers. 
    230255Indeed, time filtering is no longer required on the forcing part. 
    231 The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter 
    232 (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 
     256The influence of the Asselin filter on the forcing is explicitly removed by 
     257adding a new term in the filter (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 
    233258Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 
    234259the modified formulation becomes conservative \citep{leclair.madec_OM09}. 
    235 Second, the LF-RA becomes a truly quasi -second order scheme. 
     260Second, the LF-RA becomes a truly quasi-second order scheme. 
    236261Indeed, \autoref{eq:TD_forcing} used in combination with a careful treatment of static instability 
    237262(\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) 
    … …  
    245270even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term. 
    246271 
    247 \begin{figure}[!t] 
     272\begin{figure} 
    248273  \centering 
    249274  \includegraphics[width=0.66\textwidth]{Fig_MLF_forcing} 
    … …  
    271296\end{listing} 
    272297 
    273 The first time step of this three level scheme when starting from initial conditions is a forward step 
    274 (Euler time integration): 
     298The first time step of this three level scheme when starting from initial conditions is 
     299a forward step (Euler time integration): 
    275300\[ 
    276301  % \label{eq:TD_DOM_euler} 
    277302  x^1 = x^0 + \rdt \ \text{RHS}^0 
    278303\] 
    279 This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but 
     304This is done simply by keeping the leapfrog environment 
     305(\ie\ the \autoref{eq:TD} three level time stepping) but 
    280306setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 
    281307using half the value of a leapfrog time step ($2 \rdt$). 
    … …  
    286312running the model for $2N$ time steps in one go, 
    287313or by performing two consecutive experiments of $N$ steps with a restart. 
    288 This requires saving two time levels and many auxiliary data in the restart files in machine precision. 
     314This requires saving two time levels and many auxiliary data in 
     315the restart files in machine precision. 
    289316 
    290317Note that the time step $\rdt$, is also saved in the restart file. 
    291 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step 
    292 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting 
    293 the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model 
    294 are defined through the  \nam{run}{run} namelist variables. 
     318When restarting, if the time step has been changed, or 
     319one of the prognostic variables at \textit{before} time step is missing, 
     320an Euler time stepping scheme is imposed. 
     321A forward initial step can still be enforced by the user by 
     322setting the namelist variable \np[=0]{nn_euler}{nn\_euler}. 
     323Other options to control the time integration of the model are defined through 
     324the \nam{run}{run} namelist variables. 
     325 
    295326\gmcomment{ 
    296327add here how to force the restart to contain only one time step for operational purposes 
    … …  
    298329add also the idea of writing several restart for seasonal forecast : how is it done ? 
    299330 
    300 verify that all namelist pararmeters are truly described 
     331verify that all namelist parameters are truly described 
    301332 
    302333a word on the check of restart  ..... 
    … …  
    309340\label{subsec:TD_time} 
    310341 
    311 Options are defined through the  \nam{dom}{dom} namelist variables. 
     342Options are defined through the\nam{dom}{dom} namelist variables. 
    312343 \colorbox{yellow}{add here a few word on nit000 and nitend} 
    313344 
    314345 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} 
    315346 
    316 add a description of daymod, and the model calandar (leap-year and co) 
    317  
    318 }        %% end add 
     347add a description of daymod, and the model calendar (leap-year and co) 
     348 
     349}     %% end add 
    319350 
    320351\gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme 
Note: See TracChangeset for help on using the changeset viewer.