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 2282 for branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DOM.tex – NEMO

Ignore:
Timestamp:
2010-10-15T16:42:00+02:00 (14 years ago)
Author:
gm
Message:

ticket:#658 merge DOC of all the branches that form the v3.3 beta

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/nemo_v3_3_beta/DOC/TexFiles/Chapters/Chap_DOM.tex

    r1224 r2282  
    33% Chapter 2 Ñ Space and Time Domain (DOM) 
    44% ================================================================ 
    5 \chapter{Space and Time Domain (DOM) } 
     5\chapter{Space Domain (DOM) } 
    66\label{DOM} 
    77\minitoc 
    … …  
    1212%                  should be put outside of DOM routine (better with TRC staff and off-line 
    1313%                  tracers) 
    14 %  - daymod: definition of the time domain (nit000, nitend andd the calendar) 
    1514%  -geo2ocean:  how to switch from geographic to mesh coordinate 
    16 %  - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 
    17  
    18 \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction  
    19 which could be referred to here, would help  ==> to be added} 
    20 %%%% 
     15%     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 
    2116 
    2217 
    … …  
    2419$\ $\newline    % force a new ligne 
    2520 
    26  
    27 Having defined the continuous equations in Chap.~\ref{PE}, we need to choose a  
    28 discretization on a grid, and numerical algorithms. In the present chapter, we  
    29 provide a general description of the staggered grid used in \NEMO, and other  
    30 information relevant to the main directory routines (time stepping, main program)  
    31 as well as the DOM (DOMain) directory.  
     21Having defined the continuous equations in Chap.~\ref{PE} and chosen a time  
     22discretization Chap.~\ref{STP}, we need to choose a discretization on a grid,  
     23and numerical algorithms. In the present chapter, we provide a general description  
     24of the staggered grid used in \NEMO, and other information relevant to the main  
     25directory routines as well as the DOM (DOMain) directory.  
     26 
     27$\ $\newline    % force a new ligne 
    3228 
    3329% ================================================================ 
    … …  
    4642\begin{figure}[!tb] \label{Fig_cell}  \begin{center} 
    4743\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} 
    48 \caption{Arrangement of variables. $T$ indicates scalar points where temperature,  
     44\caption{Arrangement of variables. $t$ indicates scalar points where temperature,  
    4945salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$)  
    5046indicates vector points, and $f$ indicates vorticity points where both relative and  
    … …  
    5753Special attention has been given to the homogeneity of the solution in the three  
    5854space directions. The arrangement of variables is the same in all directions.  
    59 It consists of cells centred on scalar points ($T$, $S$, $p$, $\rho$) with vector  
     55It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector  
    6056points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}).  
    6157This is the generalisation to three dimensions of the well-known ``C'' grid in  
    … …  
    120116Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and  
    121117$k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a  
    122 variable $q$ defined at a $T$-point has its three components defined at $u$-, $v$-  
    123 and $w$-points while its Laplacien is defined at $T$-point. These operators have  
     118variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$-  
     119and $w$-points while its Laplacien is defined at $t$-point. These operators have  
    124120the following discrete forms in the curvilinear $s$-coordinate system: 
    125121\begin{equation} \label{Eq_DOM_grad} 
    126 \nabla q\equiv    \frac{1}{e_{1u} }\delta _{i+1/2} \left[ q \right]\;\,{\rm {\bf i}} 
    127          +  \frac{1}{e_{2v} }\delta _{j+1/2} \left[ q \right]\;\,{\rm {\bf j}} 
    128          +  \frac{1}{e_{3w} }\delta _{k+1/2} \left[ q \right]\;\,{\rm {\bf k}} 
     122\nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,{\rm {\bf i}} 
     123      +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,{\rm {\bf j}} 
     124      +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,{\rm {\bf k}} 
    129125\end{equation} 
    130126\begin{multline} \label{Eq_DOM_lap} 
    131 \Delta q\equiv \frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\;\left(  
    132 {\delta _i \left[ {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2}  
    133 \left[ q \right]} \right] 
    134 +\delta _j \left[ {\frac{e_{1v} e_{3v} }{e_{2v}  
    135 }\;\delta _{j+1/2} \left[ q \right]} \right]\;} \right)     \\ 
    136 +\frac{1}{e_{3T} }\delta _k \left[ {\frac{1}{e_{3w} }\;\delta _{k+1/2}  
    137 \left[ q \right]} \right] 
     127\Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }  
     128       \;\left(          \delta_i  \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] 
     129+                        \delta_j  \left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)      \\ 
     130+\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right] 
    138131\end{multline} 
    139132 
    140 Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl  
    141 components defined at $vw$-, $uw$, and $f$-points, and its divergence defined  
    142 at $T$-points: 
    143 \begin{equation} \label{Eq_DOM_curl} 
    144 \begin{split} 
    145  \nabla \times {\rm {\bf A}}\equiv \frac{1}{e_{2v} {\kern 1pt}e_{3vw}  
    146 }{\kern 1pt}\,\;\left( {\delta _{j+1/2} \left[ {e_{3w} a_3 } \right]-\delta  
    147 _{k+1/2} \left[ {e_{2v} a_2 } \right]} \right)  &\;\;{\rm {\bf i}} \\  
    148  +\frac{1}{e_{2u} {\kern 1pt}e_{3uw} }\;\left( {\delta _{k+1/2} \left[ {e_{1u} a_1 }  
    149 \right]-\delta _{i+1/2} \left[ {e_{3w} a_3 } \right]} \right)  &\;\;{\rm{\bf j}} \\  
    150  +\frac{e_{3f} }{e_{1f} {\kern 1pt}e_{2f} }\,{\kern 1pt}\;\left( {\delta  
    151 _{i+1/2} \left[ {e_{2v} a_2 } \right]-\delta _{j+12} \left[ {e_{1u} a_1 } \right]}  
    152 \right)  &\;\;{\rm {\bf k}} 
    153  \end{split} 
    154 \end{equation} 
     133Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$  
     134defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$,  
     135and $f$-points, and its divergence defined at $t$-points: 
     136\begin{eqnarray}  \label{Eq_DOM_curl} 
     137 \nabla \times {\rm {\bf A}}\equiv & 
     138      \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right)  &\ \rm{\bf i} \\  
     139 +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1  \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right)  &\ \rm{\bf j} \\  
     140 +& \frac{1}{e_{1f} \,e_{2f}    } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2  \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right)  &\ \rm{\bf k} 
     141 \end{eqnarray} 
    155142\begin{equation} \label{Eq_DOM_div} 
    156 \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} e_{3T} }\left( {\delta  
    157 _i \left[ {e_{2u} e_{3u} a_1 } \right]+\delta _j \left[ {e_{1v} e_{3v} a_2 }  
    158 \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] 
     143\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] 
     144                                                                                         +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] 
    159145\end{equation} 
    160146 
    … …  
    164150horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to:  
    165151\begin{equation*} 
    166 \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} }\left( {\delta  
    167 _i \left[ {e_{2u} a_1 } \right]+\delta _j \left[ {e_{1v}  a_2 }  
    168 \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] 
     152\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right]  
     153                                                                              +\delta_j \left[e_{1v}\, a_2 \right]  \right) 
     154                                                     +\frac{1}{e_{3t}} \delta_k \left[             a_3 \right] 
    169155\end{equation*} 
    170156 
    … …  
    172158for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): 
    173159\begin{equation} \label{DOM_bar} 
    174 \bar q   = \frac{1}{H}\int_{k^b}^{k^o} {q\;e_{3q} \,dk}  
     160\bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk}  
    175161      \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } 
    176162\end{equation} 
    … …  
    189175 
    190176It is straightforward to demonstrate that these properties are verified locally in  
    191 discrete form as soon as the scalar $q$ is taken at $T$-points and the vector  
     177discrete form as soon as the scalar $q$ is taken at $t$-points and the vector  
    192178\textbf{A} has its components defined at vector points $(u,v,w)$. 
    193179 
    … …  
    230216The array representation used in the \textsc{Fortran} code requires an integer  
    231217indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is  
    232 associated with the use of integer values for $T$-points and both integer and  
     218associated with the use of integer values for $t$-points and both integer and  
    233219integer and a half values for all the other points. Therefore a specific integer  
    234 indexing must be defined for points other than $T$-points ($i.e.$ velocity and  
     220indexing must be defined for points other than $t$-points ($i.e.$ velocity and  
    235221vorticity grid-points). Furthermore, the direction of the vertical indexing has  
    236222been changed so that the surface level is at $k=1$. 
    … …  
    242228\label{DOM_Num_Index_hor} 
    243229 
    244 The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $T$-point  
    245 and the eastward $u$-point (northward $v$-point) have the same index  
    246 (see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its  
    247 nearest northeast $f$-point have the same $i$-and $j$-indices. 
     230The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}.  
     231For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point  
     232(northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}).  
     233A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. 
    248234 
    249235% ----------------------------------- 
    … …  
    257243to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}.  
    258244The sea surface corresponds to the $w$-level $k=1$ which is the same index  
    259 as $T$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)  
     245as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)  
    260246either corresponds to the ocean floor or is inside the bathymetry while the last  
    261 $T$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that  
    262 for an increasing $k$ index, a $w$-point and the $T$-point just below have the  
     247$t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that  
     248for an increasing $k$ index, a $w$-point and the $t$-point just below have the  
    263249same $k$ index, in opposition to what is done in the horizontal plane where  
    264 it is the $T$-point and the nearest velocity points in the direction of the horizontal  
     250it is the $t$-point and the nearest velocity points in the direction of the horizontal  
    265251axis that have the same $i$ or $j$ index (compare the dashed area in  
    266252Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are  
    … …  
    309295\S\ref{LBC_mpp}). 
    310296 
     297 
     298$\ $\newline    % force a new ligne 
     299 
    311300% ================================================================ 
    312301% Domain: Horizontal Grid (mesh)  
    … …  
    341330factors in the horizontal plane as follows: 
    342331\begin{flalign*} 
    343 \lambda_T &\equiv \text{glamt}= \lambda(i)     & \varphi_T &\equiv \text{gphit} = \varphi(j)\\ 
     332\lambda_t &\equiv \text{glamt}= \lambda(i)     & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ 
    344333\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ 
    345334\lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ 
    … …  
    347336\end{flalign*} 
    348337\begin{flalign*} 
    349 e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& 
    350 e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ 
     338e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |& 
     339e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\ 
    351340e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |& 
    352341e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ 
    … …  
    359348considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with  
    360349all universal constants). Note that the horizontal position of and scale factors  
    361 at $w$-points are exactly equal to those of $T$-points, thus no specific arrays  
     350at $w$-points are exactly equal to those of $t$-points, thus no specific arrays  
    362351are defined at $w$-points.  
    363352 
    364353Note that the definition of the scale factors ($i.e.$ as the analytical first derivative  
    365354of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is  
    366 specific to the \NEMO model \citep{Marti1992}. As an example, $e_{1T}$ is defined  
    367 locally at a $T$-point, whereas many other models on a C grid choose to define  
     355specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined  
     356locally at a $t$-point, whereas many other models on a C grid choose to define  
    368357such a scale factor as the distance between the $U$-points on each side of the  
    369 $T$-point. Relying on an analytical transformation has two advantages: firstly, there  
     358$t$-point. Relying on an analytical transformation has two advantages: firstly, there  
    370359is no ambiguity in the scale factors appearing in the discrete equations, since they  
    371360are first introduced in the continuous equations; secondly, analytical transformations  
    … …  
    380369in the vertical, and (b) analytically derived grid-point position and scale factors. For  
    381370both grids here,  the same $w$-point depth has been chosen but in (a) the  
    382 $T$-points are set half way between $w$-points while in (b) they are defined from  
     371$t$-points are set half way between $w$-points while in (b) they are defined from  
    383372an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$.  
    384373Note the resulting difference between the value of the grid-size $\Delta_k$ and  
    … …  
    397386\begin{description} 
    398387\item[\jp{jphgr\_mesh}=0]  The most general curvilinear orthogonal grids. 
    399 The coordinates and their first derivatives with respect to $i$ and $j$ are  
    400 provided in a file, read in \rou{hgr\_read} subroutine of the domhgr module. 
     388The coordinates and their first derivatives with respect to $i$ and $j$ are provided 
     389in a input file (\ifile{coordinates}), read in \rou{hgr\_read} subroutine of the domhgr module. 
    401390\item[\jp{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below).  
    402391For other analytical grids, the \mdl{domhgr} module must be modified by the user.  
    … …  
    422411and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide  
    423412an approximate starting latitude: the real latitude will be recalculated analytically,  
    424 in order to ensure that the equator corresponds to line passing through $T$-  
     413in order to ensure that the equator corresponds to line passing through $t$-  
    425414and $u$-points.   
    426415 
    … …  
    431420and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively.  
    432421The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero  
    433 with the first $T$-point. The meridional coordinate (gphi. arrays) is in kilometers,  
    434 and the second $T$-point corresponds to coordinate $gphit=0$. The input  
     422with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers,  
     423and the second $t$-point corresponds to coordinate $gphit=0$. The input  
    435424parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference  
    436425latitude for computation of the Coriolis parameter. In the case of the beta plane,  
    … …  
    447436%        Grid files 
    448437% ------------------------------------------------------------------------------------------------------------- 
    449 \subsection{Grid files} 
     438\subsection{Output Grid files} 
    450439\label{DOM_hgr_files} 
    451440 
    452441All the arrays relating to a particular ocean model configuration (grid-point  
    453 position, scale factors, masks) can be saved in files if $\np{nmsh} \not= 0$  
     442position, scale factors, masks) can be saved in files if $\np{nn\_msh} \not= 0$  
    454443(namelist parameter). This can be particularly useful for plots and off-line  
    455444diagnostics. In some cases, the user may choose to make a local modification  
    … …  
    458447happens to be too wide due to insufficient model resolution). An example  
    459448is Gibraltar Strait in the ORCA2 configuration. When such modifications are done,  
    460 the output grid written when $\np{nmsh} \not=0$ is no more equal to the input grid. 
     449the output grid written when $\np{nn\_msh} \not=0$ is no more equal to the input grid. 
     450 
     451$\ $\newline    % force a new ligne 
    461452 
    462453% ================================================================ 
    … …  
    467458\label{DOM_zgr} 
    468459%-----------------------------------------nam_zgr & namdom------------------------------------------- 
    469 \namdisplay{nam_zgr}  
     460\namdisplay{namzgr}  
    470461\namdisplay{namdom}  
    471462%------------------------------------------------------------------------------------------------------------- 
    … …  
    510501Contrary to the horizontal grid, the vertical grid is computed in the code and no  
    511502provision is made for reading it from a file. The only input file is the bathymetry  
    512 (in meters)\footnote{N.B. in full step $z$-coordinate, a \textit{bathy\_level} file can  
    513 replace the \textit{bathy\_meter} file, so that the computation of the number of  
    514 wet ocean point in each water column is by-passed}. After reading the bathymetry,  
    515 the algorithm for vertical grid definition differs between the different options: 
     503(in meters) (\ifile{bathy\_meter})  
     504\footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the  
     505\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point  
     506in each water column is by-passed}.  
     507After reading the bathymetry, the algorithm for vertical grid definition differs  
     508between the different options: 
    516509\begin{description} 
    517510\item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$. 
    … …  
    535528C-Pre-Processor (CPP) key. To improve the code readability while providing this  
    536529flexibility, the vertical coordinate and scale factors are defined as functions of  
    537 $(i,j,k)$ with "fs" as prefix (examples: \textit{fsdeptht, fse3t,} etc) that can be either  
     530$(i,j,k)$ with "fs" as prefix (examples: \textit{fsdept, fse3t,} etc) that can be either  
    538531three-dimensional arrays, or a one dimensional array when \key{zco} is defined.  
    539532These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory.  
    … …  
    548541 
    549542Three options are possible for defining the bathymetry, according to the  
    550 namelist variable \np{ntopo}:  
     543namelist variable \np{nn\_bathy}:  
    551544\begin{description} 
    552 \item[\np{ntopo} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$  
     545\item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$  
    553546is given by the coordinate transformation. The domain can either be a closed  
    554547basin or a periodic channel depending on the parameter \jp{jperio}.  
    555 \item[\np{ntopo} = -1] a domain with a bump of topography one third of the  
     548\item[\np{nn\_bathy} = -1] a domain with a bump of topography one third of the  
    556549domain width at the central latitude. This is meant for the "EEL-R5" configuration,  
    557550a periodic or open boundary channel with a seamount.  
    558 \item[\np{ntopo} = 1] read a bathymetry. The bathymetry file (Netcdf format)  
     551\item[\np{nn\_bathy} = 1] read a bathymetry. The \ifile{bathy\_meter} file (Netcdf format)  
    559552provides the ocean depth (positive, in meters) at each grid point of the model grid.  
    560553The bathymetry is usually built by interpolating a standard bathymetry product  
    … …  
    563556(all levels are masked). 
    564557\end{description} 
    565  
    566 When using the rigid lid approximation (\key{dynspg\_rl} is defined) isolated land  
    567 masses (islands) must be identified by negative integers in the input bathymetry file  
    568 (see \S\ref{MISC_solisl}). 
    569558 
    570559When a global ocean is coupled to an atmospheric model it is better to represent  
    … …  
    593582 
    594583The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$  
    595 and $gdepw_0$ for $T$- and $w$-points, respectively. As indicated on  
     584and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on  
    596585Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the  
    597 ocean surface. There are at most \jp{jpk}-1 $T$-points inside the ocean, the  
    598 additional $T$-point at $jk=jpk$ is below the sea floor and is not used.  
    599 The vertical location of $w$- and $T$-levels is defined from the analytic expression  
     586ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the  
     587additional $t$-point at $jk=jpk$ is below the sea floor and is not used.  
     588The vertical location of $w$- and $t$-levels is defined from the analytic expression  
    600589of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the  
    601590vertical scale factors. The user must provide the analytical expression of both  
    … …  
    663652\begin{center} \begin{tabular}{c||r|r|r|r} 
    664653\hline 
    665 \textbf{LEVEL}& \textbf{GDEPT}& \textbf{GDEPW}& \textbf{E3T }& \textbf{E3W  } \\ \hline 
     654\textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline 
    6666551  &  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline 
    6676562  &  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline 
    … …  
    728717Two variables in the namdom namelist are used to define the partial step  
    729718vertical grid. The mimimum water thickness (in meters) allowed for a cell  
    730 partially filled with bathymetry at level jk is the minimum of \np{e3zpsmin}  
    731 (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{e3zps\_rat}$ (a fraction,  
     719partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min}  
     720(thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction,  
    732721usually 10\%, of the default thickness $e_{3t}(jk)$). 
    733722 
    … …  
    738727% ------------------------------------------------------------------------------------------------------------- 
    739728\subsection   [$s$-coordinate (\np{ln\_sco})] 
    740          {$s$-coordinate (\np{ln\_sco}=true)} 
     729           {$s$-coordinate (\np{ln\_sco}=true)} 
    741730\label{DOM_sco} 
    742731%------------------------------------------nam_zgr_sco--------------------------------------------------- 
    743 \namdisplay{nam_zgr_sco}  
     732\namdisplay{namzgr_sco}  
    744733%-------------------------------------------------------------------------------------------------------------- 
    745734In $s$-coordinate (\key{sco} is defined), the depth and thickness of the model  
    … …  
    752741\end{split} 
    753742\end{equation} 
    754 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $T$-point  
     743where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point  
    755744location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea  
    756745surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean  
    757746depth, since a mixed step-like and bottom-following representation of the  
    758747topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided  
    759 (\hf{zgr\_s} file) $h$ is a smooth envelope bathymetry and steps are used to represent  
    760 sharp bathymetric gradients. 
    761  
    762 A new flexible stretching function, modified from \citet{Song1994} is provided as an example: 
     748(\rou{zgr\_sco} routine, see \mdl{domzgr}) $h$ is a smooth envelope bathymetry  
     749and steps are used to represent sharp bathymetric gradients. 
     750 
     751A new flexible stretching function, modified from \citet{Song_Haidvogel_JCP94} is provided as an example: 
    763752\begin{equation} \label{DOM_sco_function} 
    764753\begin{split} 
    … …  
    801790Whatever the vertical coordinate used, the model offers the possibility of  
    802791representing the bottom topography with steps that follow the face of the  
    803 model cells (step like topography) \citep{Madec1996}. The distribution of  
     792model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of  
    804793the steps in the horizontal is defined in a 2D integer array, mbathy, which  
    805794gives the number of ocean levels ($i.e.$ those that are not masked) at each  
    806 $T$-point. mbathy is computed from the meter bathymetry using the definiton of  
    807 gdept as the number of $T$-points which gdept $\leq$ bathy. Note that in version  
    808 NEMO v2.3, the user still has to provide the "level" bathymetry in a NetCDF 
    809 file when using the full step option (\np{ln\_zco}), rather than the bathymetry  
    810 in meters: both will be allowed in future versions. 
     795$t$-point. mbathy is computed from the meter bathymetry using the definiton of  
     796gdept as the number of $t$-points which gdept $\leq$ bathy.  
    811797 
    812798Modifications of the model bathymetry are performed in the \textit{bat\_ctl}  
    813799routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points  
    814800that do not communicate with another ocean point at the same level are eliminated. 
    815  
    816 In the case of the rigid-lid approximation when islands occur in the computational  
    817 domain (\np{ln\_dynspg\_rl}=.true. and \key{island} is defined), the \textit{mbathy}  
    818 array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the  
    819 following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T$-points are  
    820 land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $T$-points are land points  
    821 on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points  
    822 are ocean points, the others are points below the ocean floor.  
    823  
    824 This is used to compute the island barotropic stream function used in the rigid lid  
    825 computation (see \S\ref{MISC_solisl}). 
    826801 
    827802From the \textit{mbathy} array, the mask fields are defined as follows: 
    … …  
    848823\gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  } 
    849824%%% 
    850  
    851 % ================================================================ 
    852 % Time Discretisation 
    853 % ================================================================ 
    854 \section{Time Discretisation} 
    855 \label{DOM_nxt} 
    856  
    857 The time stepping used in \NEMO is a three level scheme that can be  
    858 represented as follows: 
    859 \begin{equation} \label{Eq_DOM_nxt} 
    860    x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} 
    861 \end{equation}  
    862 where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the  
    863 corresponding time evolution equation; $\Delta t$ is the time step; and the  
    864 superscripts indicate the time at which a quantity is evaluated. Each term of the  
    865 RHS is evaluated at a specific time step(s) depending on the physics with which  
    866 it is associated.  
    867  
    868 The choice of the time step used for this evaluation is discussed below as  
    869 well as the implications in terms of starting or restarting a model  
    870 simulation. Note that the time stepping is generally performed in a one step  
    871 operation. With such a complex and nonlinear system of equations it would be  
    872 dangerous to let a prognostic variable evolve in time for each term separately. 
    873  
    874 The three level scheme requires three arrays for each prognostic variables.  
    875 For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array,  
    876 although referred to as $x_a$ (after) in the code, is usually not the variable at  
    877 the next time step; but rather it is used to store the time derivative (RHS in  
    878 \eqref{Eq_DOM_nxt}) prior to time-stepping the equation. Generally, the time  
    879 stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt}  
    880 modules, except for implicit vertical diffusion or sea surface height when  
    881 time-splitting options are used. 
    882  
    883 % ------------------------------------------------------------------------------------------------------------- 
    884 %        Non-Diffusive Part---Leapfrog Scheme 
    885 % ------------------------------------------------------------------------------------------------------------- 
    886 \subsection{Non-Diffusive Part --- Leapfrog Scheme} 
    887 \label{DOM_nxt_leap_frog} 
    888  
    889 The time stepping used for non-diffusive processes is the well-known  
    890 leapfrog scheme. It is a time centred scheme, i.e. the RHS is evaluated at  
    891 time step $t$, the now time step. It is only used for non-diffusive terms,  
    892 that is momentum and tracer advection, pressure gradient, and Coriolis  
    893 terms. This scheme is widely used for advective processes in low-viscosity  
    894 fluids. It is an efficient method that achieves second-order accuracy with  
    895 just one right hand side evaluation per time step. Moreover, it does not  
    896 artificially damp linear oscillatory motion nor does it produce instability  
    897 by amplifying the oscillations. These advantages are somewhat diminished by  
    898 the large phase-speed error of the leapfrog scheme, and the unsuitability of  
    899 leapfrog differencing for the representation of diffusive and Rayleigh  
    900 damping processes. However, the most serious problem associated with the  
    901 leapfrog scheme is a high-frequency computational noise called  
    902 "time-splitting" \citep{Haltiner1980} that develops when the method  
    903 is used to model non linear fluid dynamics: the even and odd time steps tend  
    904 to diverge into a physical and a computational mode. Time splitting can  
    905 be controlled through the use of an Asselin time filter (first designed by  
    906 \citep{Robert1966} and more comprehensively studied by \citet{Asselin1972}),  
    907 or by periodically reinitialising the leapfrog solution through a single  
    908 integration step with a two-level scheme. In \NEMO we follow the first  
    909 strategy: 
    910 \begin{equation} \label{Eq_DOM_nxt_asselin} 
    911 x_F^t  = x^t + \gamma \, \left[ x_f^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] 
    912 \end{equation}  
    913 where the subscript $f$ denotes filtered values and $\gamma$ is the Asselin  
    914 coefficient. $\gamma$ is initialized as \np{atfp} (namelist parameter).  
    915 Its default value is \np{atfp}=0.1.  This default value causes a significant dissipation  
    916 of high frequency motions. Recommended values in idealized studies of shallow  
    917 water turbulence are two orders of magnitude smaller (\citep{Farge1987}).  
    918 Both strategies do, nevertheless, degrade the accuracy of the calculation from  
    919 second to first order. The leapfrog scheme combined with a Robert-Asselin  
    920 time filter has been preferred to other time differencing schemes such as  
    921 predictor corrector or trapezoidal schemes, because the user has an explicit  
    922 and simple control of the magnitude of the time diffusion of the scheme.  
    923 In association with the 2nd order centred space discretisation of the  
    924 advective terms in the momentum and tracer equations, it avoids implicit  
    925 numerical diffusion in both the time and space discretisations of the  
    926 advective term: they are both set explicitly by the user through the Robert-Asselin  
    927 filter parameter and the viscous and diffusive coefficients. 
    928  
    929 \gmcomment{ 
    930 %gm - reflexion about leapfrog: ongoing work with Matthieu Leclair 
    931 % to be updated latter with addition of new time stepping strategy 
    932 \colorbox{yellow}{Note}:  
    933 The Robert-Asselin time filter slightly departs from a simple second order time  
    934 diffusive operator computed with a forward time stepping due to the presence of  
    935 $x_f^{t-\Delta t}$ in the right hand side of  \ref{Eq_DOM_nxt_asselin}. The original  
    936 willing of Robert1966 and Asselin1972 was to design a time filter that allow much  
    937 larger parameter than 0.5.   is due to computer saving consideration. In the original  
    938 asselin filter, $x^{t-\Delta t}$ is used instead: 
    939  \begin{equation} \label{Eq_DOM_nxt_asselin_true} 
    940 x_f^t  = x^t + \gamma \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] 
    941 \end{equation}  
    942 Applying a "true" Asselin time filter is nothing more than adding a harmonic  
    943 diffusive operator in time. Indeed, equations \ref{Eq_DOM_nxt} and  
    944 \ref{Eq_DOM_nxt_asselin_true} can be rewritten together as: 
    945 \begin{equation} \label{Eq_DOM_nxt2} 
    946 \begin{split} 
    947   \frac{ x^{t+\Delta t} - x^{t-\Delta t} } { 2 \,\Delta t }  
    948   &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t} + \frac{ x_f^t  - x^t }{2 \,\Delta t} \\ 
    949   &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t}  
    950     + \gamma\ \frac{  \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] }{2 \,\Delta t}  \\ 
    951   &=  \text{RHS}_x^{t-\Delta t,t,t+\Delta t}  
    952   + 2 \Delta t \ \gamma \ \frac{1}{{2 \Delta t}^2}   
    953    \,\delta_{t-1}\,\left[ \delta_{t+1/2}\left[ x^t \right] \right]  
    954   \end{split} 
    955 \end{equation}  
    956 expressing this again in a continuous form, one finds that the Asselin filter leads to : 
    957 \begin{equation} \label{Eq_DOM_nxt3} 
    958   \frac{ \partial x} { \partial t } =  \text{RHS} + 2 \,\Delta t \ \gamma \ \frac{ {\partial}^2 x}{ \partial t ^2 } 
    959 \end{equation}  
    960  
    961 Equations  \ref{Eq_DOM_nxt2} and \ref{Eq_DOM_nxt3} suggest several remarks.  
    962 First the Asselin filter is definitively a second order time diffusive operator which is  
    963 evaluated at centered time step. The magnitude of this diffusion is proportional to  
    964 the time step (with $\gamma$ usually taken between $10^{-1}$ to $10^{-3}$).  
    965 Second, this term has to be taken into account in all budgets of the equations  
    966 (mass, heat content, salt content, kinetic energy...). Nevertheless,we stress here  
    967 that it is small and does not create systematic biases. Indeed let us evaluate how  
    968 it contributes to the time evolution of $x$ between $t_o$ and $t_1$: 
    969 \begin{equation} \label{Eq_DOM_nxt4} 
    970 \begin{split} 
    971  t_1-t_o &= \int_{t_o}^{t_1} \frac{ \partial x} { \partial t } dt \\ 
    972       &= \int_{t_o}^{t_1} \text{RHS} dt + 2 \,\Delta t \ \gamma \left( 
    973         \left. \frac{ \partial x}{ \partial t } \right|_1  
    974       - \left. \frac{ \partial x}{ \partial t } \right|_o  \right) 
    975  \end{split} 
    976 \end{equation}  
    977 } 
    978  
    979 Alternative time stepping schemes are currently under investigation. 
    980  
    981 % ------------------------------------------------------------------------------------------------------------- 
    982 %        Diffusive Part---Forward or Backward Scheme 
    983 % ------------------------------------------------------------------------------------------------------------- 
    984 \subsection{Diffusive Part --- Forward or Backward Scheme} 
    985 \label{DOM_nxt_forward_imp} 
    986  
    987 The leapfrog differencing scheme is unsuitable for the representation of  
    988 diffusive and damping processes. For a tendancy $D_x$, representing a  
    989 diffusive term or a restoring term to a tracer climatology  
    990 (when present, see \S~\ref{TRA_dmp}), a forward time differencing scheme 
    991  is used : 
    992 \begin{equation} \label{Eq_DOM_nxt_euler} 
    993    x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ {D_x}^{t-\Delta t} 
    994 \end{equation}  
    995  
    996 This is diffusive in time and conditionally stable. For example, the  
    997 conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}: 
    998 \begin{equation} \label{Eq_DOM_nxt_euler_stability} 
    999 A^h < \left\{ 
    1000 \begin{aligned} 
    1001                     &\frac{e^2}{  8 \; \Delta t }  &&\quad \text{laplacian diffusion}  \\ 
    1002                     &\frac{e^4}{64 \; \Delta t }   &&\quad \text{bilaplacian diffusion}  
    1003             \end{aligned} 
    1004 \right. 
    1005 \end{equation} 
    1006 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is  
    1007 the mixing coefficient. The linear constraint \eqref{Eq_DOM_nxt_euler_stability}  
    1008 is a necessary condition, but not sufficient. If it is not satisfied, even mildly,  
    1009 then the model soon becomes wildly unstable. The instability can be removed  
    1010 by either reducing the length of the time steps or reducing the mixing coefficient. 
    1011  
    1012 For the vertical diffusion terms, a forward time differencing scheme can be  
    1013 used, but usually the numerical stability condition implies a strong  
    1014 constraint on the time step. Two solutions are available in \NEMO to overcome  
    1015 the stability constraint: $(a)$ a forward time differencing scheme using a  
    1016 time splitting technique (\np{ln\_zdfexp}=.true.) or $(b)$ a backward (or implicit)  
    1017 time differencing scheme by \np{ln\_zdfexp}=.false.). In $(a)$, the master  
    1018 time step $\Delta $t is cut into $N$ fractional time steps so that the  
    1019 stability criterion is reduced by a factor of $N$. The computation is done as  
    1020 follows: 
    1021 \begin{equation} \label{Eq_DOM_nxt_ts} 
    1022 \begin{split} 
    1023 & u_\ast ^{t-\Delta t} = u^{t-\Delta t}   \\ 
    1024 & u_\ast ^{t-\Delta t+L\frac{2\Delta t}{N}}=u_\ast ^{t-\Delta t+\left( {L-1}  
    1025 \right)\frac{2\Delta t}{N}}+\frac{2\Delta t}{N}\;\text{DF}^{t-\Delta t+\left( {L-1} \right)\frac{2\Delta t}{N}} 
    1026         \quad \text{for $L=1$ to $N$}      \\ 
    1027 & u^{t+\Delta t} = u_\ast^{t+\Delta t} 
    1028 \end{split} 
    1029 \end{equation} 
    1030 with DF a vertical diffusion term. The number of fractional time steps, $N$, is given  
    1031 by setting \np{n\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally  
    1032 stable but diffusive. It can be written as follows: 
    1033 \begin{equation} \label{Eq_DOM_nxt_imp} 
    1034    x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \  \text{RHS}_x^{t+\Delta t} 
    1035 \end{equation}  
    1036  
    1037 This scheme is rather time consuming since it requires a matrix inversion,  
    1038 but it becomes attractive since a splitting factor of 3 or more is needed  
    1039 for the forward time differencing scheme. For example, the finite difference  
    1040 approximation of the temperature equation is: 
    1041 \begin{equation} \label{Eq_DOM_nxt_imp_zdf} 
    1042 \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3T} }\delta  
    1043 _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]}  
    1044 \right] 
    1045 \end{equation} 
    1046 where RHS is the right hand side of the equation except for the vertical diffusion term.  
    1047 We rewrite \eqref{Eq_DOM_nxt_imp} as: 
    1048 \begin{equation} \label{Eq_DOM_nxt_imp_mat} 
    1049 -c(k+1)\;u^{t+1}(k+1)+d(k)\;u^{t+1}(k)-\;c(k)\;u^{t+1}(k-1) \equiv b(k) 
    1050 \end{equation} 
    1051 where  
    1052 \begin{align*}  
    1053  c(k) &= A_w^{vm} (k) \, / \, e_{3uw} (k)     \\ 
    1054  d(k) &= e_{3u} (k)       \, / \, (2\Delta t) + c_k + c_{k+1}    \\ 
    1055  b(k) &= e_{3u} (k) \; \left( u^{t-1}(k) \, / \, (2\Delta t) + \text{RHS} \right)     
    1056 \end{align*} 
    1057  
    1058 \eqref{Eq_DOM_nxt_imp_mat} is a linear system of equations which associated  
    1059 matrix is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal  
    1060 term is greater than the sum of the two extra-diagonal terms, therefore a special  
    1061 adaptation of the Gauss elimination procedure is used to find the solution  
    1062 (see for example \citet{Richtmyer1967}). 
    1063  
    1064 % ------------------------------------------------------------------------------------------------------------- 
    1065 %        Start/Restart strategy 
    1066 % ------------------------------------------------------------------------------------------------------------- 
    1067 \subsection{Start/Restart strategy} 
    1068 \label{DOM_nxt_rst} 
    1069 %--------------------------------------------namrun------------------------------------------- 
    1070 \namdisplay{namrun}          
    1071 %-------------------------------------------------------------------------------------------------------------- 
    1072  
    1073 The first time step of this three level scheme when starting from initial conditions  
    1074 is a forward step (Euler time integration):  
    1075 \begin{equation} \label{Eq_DOM_euler} 
    1076    x^1 = x^0 + \Delta t \ \text{RHS}^0 
    1077 \end{equation} 
    1078  
    1079 It is also possible to restart from a previous computation, by using a  
    1080 restart file. The restart strategy is designed to ensure perfect  
    1081 restartability of the code: the user should obtain the same results to  
    1082 machine precision either by running the model for $2N$ time steps in one go,  
    1083 or by performing two consecutive experiments of $N$ steps with a restart.  
    1084 This requires saving two time levels and many auxiliary data in the restart  
    1085 files in machine precision.  
    1086  
    1087 Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure  
    1088 gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be  
    1089 added in the restart file to ensure an exact restartability. This is done only optionally  
    1090 via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of  
    1091 restart file can be obtained when the restartability is not a key issue (operational  
    1092 oceanography or ensemble simulation for seasonal forcast). 
    1093 %%% 
    1094 \gmcomment{add here how to force the restart to contain only one time step for operational purposes} 
    1095 %%% 
    1096  
    1097 \gmcomment{       % add a subsection here   
    1098  
    1099 %------------------------------------------------------------------------------------------------------------- 
    1100 %        Time Domain 
    1101 % ------------------------------------------------------------------------------------------------------------- 
    1102 \subsection{Time domain} 
    1103 \label{DOM_nxt_time} 
    1104  
    1105  \colorbox{yellow}{add here a few word on nit000 and nitend} 
    1106  
    1107  \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} 
    1108  
    1109 }        %% end add 
    1110  
Note: See TracChangeset for help on using the changeset viewer.