Changeset 12178 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO/subfiles/chap_DOM.tex
- Timestamp:
- 2019-12-11T12:02:38+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc
- Property svn:ignore deleted
-
Property
svn:externals
set to
^/utils/badges badges
^/utils/logos logos
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex
- Property svn:ignore deleted
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO
-
Property
svn:externals
set to
^/utils/figures/NEMO figures
-
Property
svn:externals
set to
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO/subfiles
-
Property
svn:ignore
set to
*.ind
*.ilg
-
Property
svn:ignore
set to
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/doc/latex/NEMO/subfiles/chap_DOM.tex
r11123 r12178 2 2 3 3 \begin{document} 4 % ================================================================ 5 % Chapter 2 ——— Space and Time Domain (DOM) 6 % ================================================================ 4 7 5 \chapter{Space Domain (DOM)} 8 6 \label{chap:DOM} 9 7 10 \minitoc 11 12 % Missing things: 13 % - istate: description of the initial state ==> this has to be put elsewhere.. 14 % perhaps in MISC ? By the way the initialisation of T S and dynamics 15 % should be put outside of DOM routine (better with TRC staff and off-line 16 % tracers) 17 % -geo2ocean: how to switch from geographic to mesh coordinate 18 % - domclo: closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 19 20 \newpage 21 22 Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, 23 we need to choose a discretization on a grid, and numerical algorithms. 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 16 17 \thispagestyle{plain} 18 19 \chaptertoc 20 21 \paragraph{Changes record} ~\\ 22 23 {\footnotesize 24 \begin{tabularx}{0.8\textwidth}{l||X|X} 25 Release & 26 Author(s) & 27 Modifications \\ 28 \hline 29 {\em 4.0 } & 30 {\em Simon M\"{u}ller \& Andrew Coward \newline \newline 31 Simona Flavoni and Tim Graham } & 32 {\em Compatibility changes: many options moved to external domain configuration tools 33 (see \autoref{apdx:DOMCFG}). \newline 34 Updates } \\ 35 {\em 3.6 } & 36 {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec } & 37 {\em Updates } \\ 38 {\em $\leq$ 3.4 } & 39 {\em Gurvan Madec and S\'{e}bastien Masson } & 40 {\em First version } 41 \end{tabularx} 42 } 43 44 \clearpage 45 46 Having defined the continuous equations in \autoref{chap:MB} and 47 chosen a time discretisation \autoref{chap:TD}, 48 we need to choose a grid for spatial discretisation and related numerical algorithms. 24 49 In the present chapter, we provide a general description of the staggered grid used in \NEMO, 25 and other information relevant to the main directory routines as well as the DOM (DOMain) directory. 26 27 % ================================================================ 28 % Fundamentals of the Discretisation 29 % ================================================================ 50 and other relevant information about the DOM (DOMain) source code modules. 51 52 %% ================================================================================================= 30 53 \section{Fundamentals of the discretisation} 31 54 \label{sec:DOM_basics} 32 55 33 % ------------------------------------------------------------------------------------------------------------- 34 % Arrangement of Variables 35 % ------------------------------------------------------------------------------------------------------------- 56 %% ================================================================================================= 36 57 \subsection{Arrangement of variables} 37 58 \label{subsec:DOM_cell} 38 59 39 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 40 \begin{figure}[!tb] 41 \begin{center} 42 \includegraphics[]{Fig_cell} 43 \caption{ 44 \protect\label{fig:cell} 45 Arrangement of variables. 46 $t$ indicates scalar points where temperature, salinity, density, pressure and 47 horizontal divergence are defined. 48 $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and 49 planetary vorticities are defined. 50 } 51 \end{center} 60 \begin{figure} 61 \centering 62 \includegraphics[width=0.33\textwidth]{Fig_cell} 63 \caption[Arrangement of variables in the unit cell of space domain]{ 64 Arrangement of variables in the unit cell of space domain. 65 $t$ indicates scalar points where 66 temperature, salinity, density, pressure and horizontal divergence are defined. 67 $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where 68 both relative and planetary vorticities are defined.} 69 \label{fig:DOM_cell} 52 70 \end{figure} 53 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 54 55 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, 56 centred second-order finite difference approximation. 57 Special attention has been given to the homogeneity of the solution in the three space directions. 71 72 The numerical techniques used to solve the Primitive Equations in this model are based on 73 the traditional, centred second-order finite difference approximation. 74 Special attention has been given to the homogeneity of the solution in the three spatial directions. 58 75 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 points $(u, v, w)$ defined in 60 the centre of each face of the cells (\autoref{fig:cell}). 61 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 62 \citep{mesinger.arakawa_bk76}. 63 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 64 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 65 66 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 67 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 68 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. 69 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 70 the grid-point where the scale factors are defined. 71 Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. 76 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with 77 vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). 78 This is the generalisation to three dimensions of the well-known ``C'' grid in 79 Arakawa's classification \citep{mesinger.arakawa_bk76}. 80 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each 81 vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying 82 the $\zeta$ and $f$-points. 83 84 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by 85 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 86 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on 87 \autoref{tab:DOM_cell}. 88 In all the following, 89 subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where 90 the scale factors are defined. 91 Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 72 92 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 73 $\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity. 74 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 75 while the scale factors are chosen equal to their local analytical value. 93 $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 94 Discrete partial derivatives are formulated by 95 the traditional, centred second order finite difference approximation while 96 the scale factors are chosen equal to their local analytical value. 76 97 An important point here is that the partial derivative of the scale factors must be evaluated by 77 98 centred finite difference approximation, not from their analytical expression. 78 This preserves the symmetry of the discrete set of equations and therefore satisfies many of79 the continuous properties (see \autoref{apdx:C}).99 This preserves the symmetry of the discrete set of equations and 100 therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). 80 101 A similar, related remark can be made about the domain size: 81 when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors 82 (see \autoref{eq:DOM_bar} in the next section). 83 84 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 85 \begin{table}[!tb] 86 \begin{center} 87 \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} 88 \hline 89 T & $i $ & $j $ & $k $ \\ 90 \hline 91 u & $i + 1/2$ & $j $ & $k $ \\ 92 \hline 93 v & $i $ & $j + 1/2$ & $k $ \\ 94 \hline 95 w & $i $ & $j $ & $k + 1/2$ \\ 96 \hline 97 f & $i + 1/2$ & $j + 1/2$ & $k $ \\ 98 \hline 99 uw & $i + 1/2$ & $j $ & $k + 1/2$ \\ 100 \hline 101 vw & $i $ & $j + 1/2$ & $k + 1/2$ \\ 102 \hline 103 fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 104 \hline 105 \end{tabular} 106 \caption{ 107 \protect\label{tab:cell} 108 Location of grid-points as a function of integer or integer and a half value of the column, line or level. 109 This indexing is only used for the writing of the semi -discrete equation. 110 In the code, the indexing uses integer values only and has a reverse direction in the vertical 111 (see \autoref{subsec:DOM_Num_Index}) 112 } 113 \end{center} 102 when needed, an area, volume, or the total ocean depth must be evaluated as 103 the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). 104 105 \begin{table} 106 \centering 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$ \\ 124 \hline 125 \end{tabular} 126 \caption[Location of grid-points]{ 127 Location of grid-points as a function of integer or 128 integer and a half value of the column, line or level. 129 This indexing is only used for the writing of the semi-discrete equations. 130 In the code, the indexing uses integer values only and 131 is positive downwards in the vertical with $k=1$ at the surface. 132 (see \autoref{subsec:DOM_Num_Index})} 133 \label{tab:DOM_cell} 114 134 \end{table} 115 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 116 117 % ------------------------------------------------------------------------------------------------------------- 118 % Vector Invariant Formulation 119 % ------------------------------------------------------------------------------------------------------------- 135 136 Note that the definition of the scale factors 137 (\ie\ as the analytical first derivative of the transformation that 138 results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 139 is specific to the \NEMO\ model \citep{marti.madec.ea_JGR92}. 140 As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, 141 whereas many other models on a C grid choose to define such a scale factor as 142 the distance between the $u$-points on each side of the $t$-point. 143 Relying on an analytical transformation has two advantages: 144 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 145 since they are first introduced in the continuous equations; 146 secondly, analytical transformations encourage good practice by 147 the 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}. 150 An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 151 \begin{figure} 152 \centering 153 \includegraphics[width=0.5\textwidth]{Fig_zgr_e3} 154 \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ 155 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 156 and (b) analytically derived grid-point position and scale factors. 157 For both grids here, the same $w$-point depth has been chosen but 158 in (a) the $t$-points are set half way between $w$-points while 159 in (b) they are defined from an analytical function: 160 $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 161 Note the resulting difference between the value of the grid-size $\Delta_k$ and 162 those of the scale factor $e_k$.} 163 \label{fig:DOM_zgr_e3} 164 \end{figure} 165 166 %% ================================================================================================= 120 167 \subsection{Discrete operators} 121 168 \label{subsec:DOM_operators} 122 169 123 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at124 the midpoint between them are:170 Given the values of a variable $q$ at adjacent points, 171 the differencing and averaging operators at the midpoint between them are: 125 172 \begin{alignat*}{2} 126 % \label{eq: di_mi}173 % \label{eq:DOM_di_mi} 127 174 \delta_i [q] &= & &q (i + 1/2) - q (i - 1/2) \\ 128 175 \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 … … 130 177 131 178 Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 132 Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at 133 a $t$-point has its three components defined at $u$-, $v$- and $w$-points while 134 its Laplacian is defined at $t$-point. 179 Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, 180 the gradient of a variable $q$ defined at a $t$-point has 181 its three components defined at $u$-, $v$- and $w$-points while 182 its Laplacian is defined at the $t$-point. 135 183 These operators have the following discrete forms in the curvilinear $s$-coordinates system: 136 \ [184 \begin{gather*} 137 185 % \label{eq:DOM_grad} 138 186 \nabla q \equiv \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i 139 187 + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j 140 + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 141 \] 142 \begin{multline*} 188 + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\ 143 189 % \label{eq:DOM_lap} 144 190 \Delta q \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 145 191 \; \lt[ \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) 146 + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\192 + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] 147 193 + \frac{1}{e_{3t}} 148 194 \delta_k \lt[ \frac{1 }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 149 \end{multline*} 150 151 Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 152 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 195 \end{gather*} 196 197 Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, 198 a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has 199 its three curl components defined at $vw$-, $uw$, and $f$-points, and 153 200 its divergence defined at $t$-points: 154 \begin{multline }201 \begin{multline*} 155 202 % \label{eq:DOM_curl} 156 203 \nabla \times \vect A \equiv \frac{1}{e_{2v} \, e_{3vw}} … … 163 210 \Big[ \delta_{i + 1/2} (e_{2v} \, a_2) 164 211 - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 165 \end{multline }166 \ begin{equation}212 \end{multline*} 213 \[ 167 214 % \label{eq:DOM_div} 168 215 \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 169 216 \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 170 217 + \frac{1}{e_{3t}} \delta_k (a_3) 171 \ end{equation}172 173 The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which174 is a masked field (i.e. equal to zero inside solid area):218 \] 219 220 The vertical average over the whole water column is denoted by an overbar and 221 is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 175 222 \begin{equation} 176 223 \label{eq:DOM_bar} … … 178 225 \end{equation} 179 226 where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 180 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over 181 all grid points of the same type in the direction indicated by the subscript (here $k$). 227 $k^b$ and $k^o$ are the bottom and surface $k$-indices, 228 and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in 229 the direction indicated by the subscript (here $k$). 182 230 183 231 In continuous form, the following properties are satisfied: … … 189 237 \end{gather} 190 238 191 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as192 the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at239 It is straightforward to demonstrate that these properties are verified locally in discrete form as 240 soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 193 241 vector points $(u,v,w)$. 194 242 195 Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. 196 Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 197 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 198 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 199 \begin{alignat}{4} 243 Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 244 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and 245 $\delta_k$) are skew-symmetric linear operators, 246 and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and 247 $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 248 \begin{alignat}{5} 200 249 \label{eq:DOM_di_adj} 201 250 &\sum \limits_i a_i \; \delta_i [b] &\equiv &- &&\sum \limits_i \delta _{ i + 1/2} [a] &b_{i + 1/2} \\ … … 204 253 \end{alignat} 205 254 206 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 255 In other words, 256 the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 207 257 $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 208 These two properties will be used extensively in the \autoref{apdx: C} to258 These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to 209 259 demonstrate integral conservative properties of the discrete formulation chosen. 210 260 211 % ------------------------------------------------------------------------------------------------------------- 212 % Numerical Indexing 213 % ------------------------------------------------------------------------------------------------------------- 261 %% ================================================================================================= 214 262 \subsection{Numerical indexing} 215 263 \label{subsec:DOM_Num_Index} 216 264 217 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 218 \begin{figure}[!tb] 219 \begin{center} 220 \includegraphics[]{Fig_index_hor} 221 \caption{ 222 \protect\label{fig:index_hor} 223 Horizontal integer indexing used in the \fortran code. 224 The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices 225 } 226 \end{center} 265 \begin{figure} 266 \centering 267 \includegraphics[width=0.33\textwidth]{Fig_index_hor} 268 \caption[Horizontal integer indexing]{ 269 Horizontal integer indexing used in the \fortran\ code. 270 The dashed area indicates the cell in which 271 variables contained in arrays have the same $i$- and $j$-indices} 272 \label{fig:DOM_index_hor} 227 273 \end{figure} 228 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 229 230 The array representation used in the \fortran code requires an integer indexing while 231 the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 232 integer values for $t$-points and both integer and integer and a half values for all the other points. 233 Therefore a specific integer indexing must be defined for points other than $t$-points 234 (\ie velocity and vorticity grid-points). 235 Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$. 236 237 % ----------------------------------- 238 % Horizontal Indexing 239 % ----------------------------------- 274 275 The array representation used in the \fortran\ code requires an integer indexing. 276 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with 277 the use of integer values for $t$-points only while 278 all the other points involve integer and a half values. 279 Therefore, a specific integer indexing has been defined for points other than $t$-points 280 (\ie\ velocity and vorticity grid-points). 281 Furthermore, the direction of the vertical indexing has been reversed and 282 the surface level set at $k = 1$. 283 284 %% ================================================================================================= 240 285 \subsubsection{Horizontal indexing} 241 286 \label{subsec:DOM_Num_Index_hor} 242 287 243 The indexing in the horizontal plane has been chosen as shown in \autoref{fig: index_hor}.288 The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}. 244 289 For an increasing $i$ index ($j$ index), 245 290 the $t$-point and the eastward $u$-point (northward $v$-point) have the same index 246 (see the dashed area in \autoref{fig:index_hor}). 247 A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. 248 249 % ----------------------------------- 250 % Vertical indexing 251 % ----------------------------------- 291 (see the dashed area in \autoref{fig:DOM_index_hor}). 292 A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices. 293 294 %% ================================================================================================= 252 295 \subsubsection{Vertical indexing} 253 296 \label{subsec:DOM_Num_Index_vertical} 254 297 255 In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in 256 the \fortran code compared to the indexing used in the semi -discrete equations and 257 given in \autoref{subsec:DOM_cell}. 258 The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below 259 (\autoref{fig:index_vert}). 260 The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while 261 the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}). 262 Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, 263 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 axis that 265 have the same $i$ or $j$ index 266 (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). 298 In the vertical, the chosen indexing requires special attention since 299 the direction of the $k$-axis in the \fortran\ code is the reverse of 300 that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 301 The sea surface corresponds to the $w$-level $k = 1$, 302 which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). 303 The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 304 the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 305 Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 306 (\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 307 in contrast to the indexing on the horizontal plane where 308 the $t$-point has the same index as the nearest velocity points in 309 the positive direction of the respective horizontal axis index 310 (compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 267 311 Since the scale factors are chosen to be strictly positive, 268 a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of 269 the discrete equations given in this documentation. 270 271 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 272 \begin{figure}[!pt] 273 \begin{center} 274 \includegraphics[]{Fig_index_vert} 275 \caption{ 276 \protect\label{fig:index_vert} 277 Vertical integer indexing used in the \fortran code. 278 Note that the $k$-axis is orientated downward. 279 The dashed area indicates the cell in which variables contained in arrays have the same $k$-index. 280 } 281 \end{center} 312 a \textit{minus sign} is included in the \fortran\ implementations of 313 \textit{all the vertical derivatives} of the discrete equations given in this manual in order to 314 accommodate the opposing vertical index directions in implementation and documentation. 315 316 \begin{figure} 317 \centering 318 \includegraphics[width=0.33\textwidth]{Fig_index_vert} 319 \caption[Vertical integer indexing]{ 320 Vertical integer indexing used in the \fortran\ code. 321 Note that the $k$-axis is oriented downward. 322 The dashed area indicates the cell in which 323 variables contained in arrays have a common $k$-index.} 324 \label{fig:DOM_index_vert} 282 325 \end{figure} 283 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 284 285 % ----------------------------------- 286 % Domain Size 287 % ----------------------------------- 288 \subsubsection{Domain size} 326 327 %% ================================================================================================= 328 \section{Spatial domain configuration} 329 \label{subsec:DOM_config} 330 331 Two typical methods are available to specify the spatial domain configuration; 332 they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in 333 namelist \nam{cfg}{cfg}. 334 335 If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, 336 the domain-specific parameters and fields are read from a NetCDF input file, 337 whose name (without its .nc suffix) can be specified as 338 the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 339 340 If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, 341 the domain-specific parameters and fields can be provided (\eg\ analytically computed) by 342 subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 343 These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, 344 and default versions that configure the spatial domain for the GYRE reference configuration are 345 present in the \path{./src/OCE/USR} directory. 346 347 In version 4.0 there are no longer any options for reading complex bathymetries and 348 performing a vertical discretisation at run-time. 349 Whilst it is occasionally convenient to have a common bathymetry file and, for example, 350 to run similar models with and without partial bottom boxes and/or sigma-coordinates, 351 supporting such choices leads to overly complex code. 352 Worse still is the difficulty of ensuring the model configurations intended to be identical are 353 indeed so when the model domain itself can be altered by runtime selections. 354 The code previously used to perform vertical discretisation has been incorporated into 355 an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 356 357 The next subsections summarise the parameter and fields related to 358 the configuration of the whole model domain. 359 These represent the minimum information that must be provided either via 360 the \np{cn_domcfg}{cn\_domcfg} file or 361 set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 362 The requirements are presented in three sections: 363 the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), 364 and the vertical grid (\autoref{subsec:DOM_zgr}). 365 366 %% ================================================================================================= 367 \subsection{Domain size} 289 368 \label{subsec:DOM_size} 290 369 291 The total size of the computational domain is set by the parameters \np{jpiglo}, 292 \np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively. 293 Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when 294 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined, 295 see \autoref{sec:LBC_mpp}). 296 297 % ================================================================ 298 % Domain: List of fields needed 299 % ================================================================ 300 \section{Needed fields} 301 \label{sec:DOM_fields} 302 The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that 303 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 304 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 305 The associated scale factors are defined using the analytical first derivative of the transformation 306 \autoref{eq:scale_factors}. 307 Necessary fields for configuration definition are: 370 The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and 371 \jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. 372 Note, that the variables \texttt{jpi} and \texttt{jpj} refer to 373 the size of each processor subdomain when the code is run in parallel using domain decomposition 374 (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). 375 376 The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, 377 and the nominal resolution through parameter \np{nn_cfg}{nn\_cfg} 378 (unless in the input file both of variables \texttt{ORCA} and \texttt{ORCA\_index} are present, 379 in which case \np{cn_cfg}{cn\_cfg} and \np{nn_cfg}{nn\_cfg} are set from these values accordingly). 380 381 The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. 382 See \autoref{sec:LBC_jperio} for details on the available options and 383 the corresponding values for \jp{jperio}. 384 385 %% ================================================================================================= 386 \subsection[Horizontal grid mesh (\textit{domhgr.F90}]{Horizontal grid mesh (\protect\mdl{domhgr})} 387 \label{subsec:DOM_hgr} 388 389 %% ================================================================================================= 390 \subsubsection{Required fields} 391 \label{sec:DOM_hgr_fields} 392 393 The explicit specification of a range of mesh-related fields are required for 394 the definition of a configuration. 395 These include: 396 397 \begin{clines} 398 int jpiglo, jpjglo, jpkglo /* global domain sizes */ 399 int jperio /* lateral global domain b.c. */ 400 double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ 401 double gphit, gphiu, gphiv, gphif /* geographic latitude */ 402 double e1t, e1u, e1v, e1f /* horizontal scale factors */ 403 double e2t, e2u, e2v, e2f /* horizontal scale factors */ 404 \end{clines} 405 406 The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to 407 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, 408 evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position. 409 The calculation of the values of the horizontal scale factor arrays in general additionally involves 410 partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, 411 evaluated for the same arguments as $\lambda$ and $\varphi$. 412 413 %% ================================================================================================= 414 \subsubsection{Optional fields} 415 416 \begin{clines} 417 /* Optional: */ 418 int ORCA, ORCA_index /* configuration name, configuration resolution */ 419 double e1e2u, e1e2v /* U and V surfaces (if grid size reduction in some straits) */ 420 double ff_f, ff_t /* Coriolis parameter (if not on the sphere) */ 421 \end{clines} 422 423 \NEMO\ can support the local reduction of key strait widths by 424 altering individual values of e2u or e1v at the appropriate locations. 425 This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points 426 (see \autoref{sec:MISC_strait} for illustrated examples). 427 The key is to reduce the faces of $T$-cell 428 (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 429 not the volume of the cells. 430 Doing otherwise can lead to numerical instability issues. 431 In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but 432 in cases where a gridsize reduction is required, 433 the unaltered surface areas at $u$ and $v$ grid points 434 (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. 435 If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and 436 the internal computation is suppressed. 437 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should 438 set the surface-area computation flag: 439 \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 440 441 \smallskip 442 Similar logic applies to the other optional fields: 443 \texttt{ff\_f} and \texttt{ff\_t} which can be used to 444 provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. 445 If present these fields will be read and used and 446 the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 447 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should 448 set the Coriolis computation flag: 449 \texttt{iff} to a non-zero value to suppress their re-computation. 450 451 Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to 452 those of $t$ points, thus no specific arrays are defined at $w$ points. 453 454 %% ================================================================================================= 455 \subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} 456 \label{subsec:DOM_zgr} 457 458 \begin{listing} 459 \nlst{namdom} 460 \caption{\forcode{&namdom}} 461 \label{lst:namdom} 462 \end{listing} 463 464 In the vertical, the model mesh is determined by four things: 465 \begin{enumerate} 466 \item the bathymetry given in meters; 467 \item the number of levels of the model (\jp{jpk}); 468 \item the analytical transformation $z(i,j,k)$ and the vertical scale factors 469 (derivatives of the transformation); and 470 \item the masking system, 471 \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. 472 \end{enumerate} 473 474 \begin{figure} 475 \centering 476 \includegraphics[width=0.5\textwidth]{Fig_z_zps_s_sps} 477 \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ 478 The ocean bottom as seen by the model: 479 \begin{enumerate*}[label=(\textit{\alph*})] 480 \item $z$-coordinate with full step, 481 \item $z$-coordinate with partial step, 482 \item $s$-coordinate: terrain following representation, 483 \item hybrid $s-z$ coordinate, 484 \item hybrid $s-z$ coordinate with partial step, and 485 \item same as (e) but in the non-linear free surface 486 (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 487 \end{enumerate*} 488 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 489 \label{fig:DOM_z_zps_s_sps} 490 \end{figure} 491 492 The choice of a vertical coordinate is made when setting up the configuration; 493 it is not intended to be an option which can be changed in the middle of an experiment. 494 The one exception to this statement being the choice of linear or non-linear free surface. 495 In v4.0 the linear free surface option is implemented as 496 a special case of the non-linear free surface. 497 This is computationally wasteful since it uses the structures for time-varying 3D metrics 498 for fields that (in the linear free surface case) are fixed. 499 However, the linear free-surface is rarely used and 500 implementing it this way means a single configuration file can support both options. 501 502 By default a non-linear free surface is used 503 (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 504 the coordinate follow the time-variation of the free surface so that 505 the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 506 When a linear free surface is assumed 507 (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 508 the vertical coordinates are fixed in time, but 509 the seawater can move up and down across the $z_0$ surface 510 (in other words, the top of the ocean in not a rigid lid). 511 512 Note that settings: 513 \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} 514 mentioned in the following sections appear to be namelist options but 515 they are no longer truly namelist options for \NEMO. 516 Their value is written to and read from the domain configuration file and 517 they should be treated as fixed parameters for a particular configuration. 518 They are namelist options for the \texttt{DOMAINcfg} tool that can be used to 519 build the configuration file and serve both to provide a record of the choices made whilst 520 building the configuration and to trigger appropriate code blocks within \NEMO. 521 These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. 522 523 \medskip 524 The decision on these choices must be made when the \np{cn_domcfg}{cn\_domcfg} file is constructed. 525 Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c): 308 526 309 527 \begin{itemize} 310 \item 311 Geographic position: 312 longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and 313 latitude with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif} 314 (all respectively at T, U, V and F point) 315 \item 316 Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t} 317 (at T and F point) 318 \item 319 Scale factors: 320 \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction), 321 \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and 322 \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\ 323 \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits), 324 \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed. 528 \item $z$-coordinate with full step bathymetry (\np[=.true.]{ln_zco}{ln\_zco}), 529 \item $z$-coordinate with partial step ($zps$) bathymetry (\np[=.true.]{ln_zps}{ln\_zps}), 530 \item Generalized, $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}). 325 531 \end{itemize} 326 327 These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in 328 \ngn{namcfg}. 329 330 \nlst{namcfg} 331 332 Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration. 333 For Reference Configurations of NEMO input domain files are supplied by NEMO System Team. 334 For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 335 They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory, 336 they provide the horizontal and vertical mesh. 337 % ------------------------------------------------------------------------------------------------------------- 338 % Needed fields 339 % ------------------------------------------------------------------------------------------------------------- 340 %\subsection{List of needed fields to build DOMAIN} 341 %\label{subsec:DOM_fields_list} 342 343 344 % ================================================================ 345 % Domain: Horizontal Grid (mesh) 346 % ================================================================ 347 \section{Horizontal grid mesh (\protect\mdl{domhgr})} 348 \label{sec:DOM_hgr} 349 350 % ------------------------------------------------------------------------------------------------------------- 351 % Coordinates and scale factors 352 % ------------------------------------------------------------------------------------------------------------- 353 \subsection{Coordinates and scale factors} 354 \label{subsec:DOM_hgr_coord_e} 355 356 The ocean mesh (\ie the position of all the scalar and vector points) is defined by 357 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 358 The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}. 359 The associated scale factors are defined using the analytical first derivative of the transformation 360 \autoref{eq:scale_factors}. 361 These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, 362 which provide the horizontal and vertical meshes, respectively. 363 This section deals with the horizontal mesh parameters. 364 365 In a horizontal plane, the location of all the model grid points is defined from 366 the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$. 367 The horizontal scale factors are calculated using \autoref{eq:scale_factors}. 368 For example, when the longitude and latitude are function of a single value 369 ($i$ and $j$, respectively) (geographical configuration of the mesh), 370 the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$, 371 and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module. 372 The model computes the grid-point positions and scale factors in the horizontal plane as follows: 532 533 Additionally, hybrid combinations of the three main coordinates are available: 534 $s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e). 535 536 A further choice related to vertical coordinate concerns 537 the presence (or not) of ocean cavities beneath ice shelves within the model domain. 538 A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that 539 the domain contains ocean cavities, 540 otherwise the top, wet layer of the ocean will always be at the ocean surface. 541 This option is currently only available for $z$- or $zps$-coordinates. 542 In the latter case, partial steps are also applied at the ocean/ice shelf interface. 543 544 Within the model, 545 the arrays describing the grid point depths and vertical scale factors are 546 three set of three dimensional arrays $(i,j,k)$ defined at 547 \textit{before}, \textit{now} and \textit{after} time step. 548 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 549 They are updated at each model time step. 550 The initial fixed reference coordinate system is held in variable names with a $\_0$ suffix. 551 When the linear free surface option is used (\np[=.true.]{ln_linssh}{ln\_linssh}), 552 \textit{before}, \textit{now} and \textit{after} arrays are initially set to 553 their reference counterpart and remain fixed. 554 555 %% ================================================================================================= 556 \subsubsection{Required fields} 557 \label{sec:DOM_zgr_fields} 558 559 The explicit specification of a range of fields related to the vertical grid are required for 560 the definition of a configuration. 561 These include: 562 563 \begin{clines} 564 int ln_zco, ln_zps, ln_sco /* flags for z-coord, z-coord with partial steps and s-coord */ 565 int ln_isfcav /* flag for ice shelf cavities */ 566 double e3t_1d, e3w_1d /* reference vertical scale factors at T and W points */ 567 double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ 568 double e3uw_0, e3vw_0 /* vertical scale factors 3D coordinate at UW and VW points */ 569 int bottom_level, top_level /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ 570 /* For reference: */ 571 float bathy_metry /* bathymetry used in setting top and bottom levels */ 572 \end{clines} 573 574 This set of vertical metrics is sufficient to describe the initial depth and thickness of 575 every gridcell in the model regardless of the choice of vertical coordinate. 576 With constant z-levels, e3 metrics will be uniform across each horizontal level. 577 In the partial step case each e3 at the \jp{bottom\_level} 578 (and, possibly, \jp{top\_level} if ice cavities are present) 579 may vary from its horizontal neighbours. 580 And, in s-coordinates, variations can occur throughout the water column. 581 With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that 582 variations occur throughout the water column with displacements related to the sea surface height. 583 These variations are typically much smaller than those arising from bottom fitted coordinates. 584 The values for vertical metrics supplied in the domain configuration file can be considered as 585 those arising from a flat sea surface with zero elevation. 586 587 The \jp{bottom\_level} and \jp{top\_level} 2D arrays define 588 the \jp{bottom\_level} and top wet levels in each grid column. 589 Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). 590 With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. 591 592 %% ================================================================================================= 593 \subsubsection{Level bathymetry and mask} 594 \label{subsec:DOM_msk} 595 596 From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: 373 597 \begin{align*} 374 \lambda_t &\equiv \text{glamt} = \lambda (i ) 375 &\varphi_t &\equiv \text{gphit} = \varphi (j ) \\ 376 \lambda_u &\equiv \text{glamu} = \lambda (i + 1/2) 377 &\varphi_u &\equiv \text{gphiu} = \varphi (j ) \\ 378 \lambda_v &\equiv \text{glamv} = \lambda (i ) 379 &\varphi_v &\equiv \text{gphiv} = \varphi (j + 1/2) \\ 380 \lambda_f &\equiv \text{glamf} = \lambda (i + 1/2) 381 &\varphi_f &\equiv \text{gphif} = \varphi (j + 1/2) \\ 382 e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j ) | 383 &e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 384 e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j ) | 385 &e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j ) | \\ 386 e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i ) \; \cos\varphi(j + 1/2) | 387 &e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | \\ 388 e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) | 389 &e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j + 1/2) | 598 tmask(i,j,k) &= 599 \begin{cases} 600 0 &\text{if $ k < top\_level(i,j)$} \\ 601 1 &\text{if $ bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 602 0 &\text{if $k > bottom\_level(i,j) $} 603 \end{cases} \\ 604 umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j, k) \\ 605 vmask(i,j,k) &= tmask(i,j,k) * tmask(i ,j + 1,k) \\ 606 fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j, k) * tmask(i,j,k) * tmask(i + 1,j, k) \\ 607 wmask(i,j,k) &= tmask(i,j,k) * tmask(i ,j,k - 1) \\ 608 \text{with~} wmask(i,j,1) &= tmask(i,j,1) 390 609 \end{align*} 391 where the last letter of each computational name indicates the grid point considered and 392 $r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants). 393 Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points, 394 thus no specific arrays are defined at $w$-points. 395 396 Note that the definition of the scale factors 397 (\ie as the analytical first derivative of the transformation that 398 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) 399 is specific to the \NEMO model \citep{marti.madec.ea_JGR92}. 400 As an example, $e_{1t}$ is defined locally at a $t$-point, 401 whereas many other models on a C grid choose to define such a scale factor as 402 the distance between the $U$-points on each side of the $t$-point. 403 Relying on an analytical transformation has two advantages: 404 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 405 since they are first introduced in the continuous equations; 406 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 407 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 408 An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. 409 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 410 \begin{figure}[!t] 411 \begin{center} 412 \includegraphics[]{Fig_zgr_e3} 413 \caption{ 414 \protect\label{fig:zgr_e3} 415 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, 416 and (b) analytically derived grid-point position and scale factors. 417 For both grids here, the same $w$-point depth has been chosen but 418 in (a) the $t$-points are set half way between $w$-points while 419 in (b) they are defined from an analytical function: 420 $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. 421 Note the resulting difference between the value of the grid-size $\Delta_k$ and 422 those of the scale factor $e_k$. 423 } 424 \end{center} 425 \end{figure} 426 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 427 428 % ------------------------------------------------------------------------------------------------------------- 429 % Choice of horizontal grid 430 % ------------------------------------------------------------------------------------------------------------- 431 \subsection{Choice of horizontal grid} 432 \label{subsec:DOM_hgr_msh_choice} 433 434 % ------------------------------------------------------------------------------------------------------------- 435 % Grid files 436 % ------------------------------------------------------------------------------------------------------------- 437 \subsection{Output grid files} 438 \label{subsec:DOM_hgr_files} 439 440 All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks) 441 can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}). 442 This can be particularly useful for plots and off-line diagnostics. 443 In some cases, the user may choose to make a local modification of a scale factor in the code. 444 This is the case in global configurations when restricting the width of a specific strait 445 (usually a one-grid-point strait that happens to be too wide due to insufficient model resolution). 446 An example is Gibraltar Strait in the ORCA2 configuration. 447 When such modifications are done, 448 the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid. 449 450 % ================================================================ 451 % Domain: Vertical Grid (domzgr) 452 % ================================================================ 453 \section{Vertical grid (\protect\mdl{domzgr})} 454 \label{sec:DOM_zgr} 455 %-----------------------------------------nam_zgr & namdom------------------------------------------- 456 % 457 %\nlst{namzgr} 458 459 \nlst{namdom} 460 %------------------------------------------------------------------------------------------------------------- 461 462 Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists. 463 In the vertical, the model mesh is determined by four things: 464 (1) the bathymetry given in meters; 465 (2) the number of levels of the model (\jp{jpk}); 466 (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 467 (4) the masking system, \ie the number of wet model levels at each 468 $(i,j)$ column of points. 469 470 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 471 \begin{figure}[!tb] 472 \begin{center} 473 \includegraphics[]{Fig_z_zps_s_sps} 474 \caption{ 475 \protect\label{fig:z_zps_s_sps} 476 The ocean bottom as seen by the model: 477 (a) $z$-coordinate with full step, 478 (b) $z$-coordinate with partial step, 479 (c) $s$-coordinate: terrain following representation, 480 (d) hybrid $s-z$ coordinate, 481 (e) hybrid $s-z$ coordinate with partial step, and 482 (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}~\forcode{= .false.}). 483 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e). 484 } 485 \end{center} 486 \end{figure} 487 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 488 489 The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, 490 must be done once of all at the beginning of an experiment. 491 It is not intended as an option which can be enabled or disabled in the middle of an experiment. 492 Three main choices are offered (\autoref{fig:z_zps_s_sps}): 493 $z$-coordinate with full step bathymetry (\np{ln\_zco}~\forcode{= .true.}), 494 $z$-coordinate with partial step bathymetry (\np{ln\_zps}~\forcode{= .true.}), 495 or generalized, $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}). 496 Hybridation of the three main coordinates are available: 497 $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}). 498 By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that 499 the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}). 500 When a linear free surface is assumed (\np{ln\_linssh}~\forcode{= .true.}), 501 the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface 502 (in other words, the top of the ocean in not a rigid-lid). 503 The last choice in terms of vertical coordinate concerns the presence (or not) in 504 the model domain of ocean cavities beneath ice shelves. 505 Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in. 506 This option is currently only available in $z$- or $zps$-coordinate, 507 and partial step are also applied at the ocean/ice shelf interface. 508 509 Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for 510 reading it from a file. 511 The only input file is the bathymetry (in meters) (\ifile{bathy\_meter}) 512 \footnote{ 513 N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file, 514 so that the computation of the number of wet ocean point in each water column is by-passed}. 515 If \np{ln\_isfcav}~\forcode{= .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing 516 the ice shelf draft (in meters) is needed. 517 518 After reading the bathymetry, the algorithm for vertical grid definition differs between the different options: 519 \begin{description} 520 \item[\textit{zco}] 521 set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$. 522 \item[\textit{zps}] 523 set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at 524 each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays. 525 \item[\textit{sco}] 526 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 527 set the three-dimensional transformation. 528 \item[\textit{s-z} and \textit{s-zps}] 529 smooth the bathymetry to fulfill the hydrostatic consistency criteria and 530 set the three-dimensional transformation $z(i,j,k)$, 531 and possibly introduce masking of extra land points to better fit the original bathymetry file. 532 \end{description} 533 %%% 534 \gmcomment{ add the description of the smoothing: envelop topography...} 535 %%% 536 537 Unless a linear free surface is used (\np{ln\_linssh}~\forcode{= .false.}), 538 the arrays describing the grid point depths and vertical scale factors are three set of 539 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 540 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 541 They are updated at each model time step using a fixed reference coordinate system which 542 computer names have a $\_0$ suffix. 543 When the linear free surface option is used (\np{ln\_linssh}~\forcode{= .true.}), \textit{before}, 544 \textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart. 545 546 % ------------------------------------------------------------------------------------------------------------- 547 % Meter Bathymetry 548 % ------------------------------------------------------------------------------------------------------------- 549 \subsection{Meter bathymetry} 550 \label{subsec:DOM_bathy} 551 552 Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy} 553 (found in \ngn{namdom} namelist): 554 \begin{description} 555 \item[\np{nn\_bathy}~\forcode{= 0}]: 556 a flat-bottom domain is defined. 557 The total depth $z_w (jpk)$ is given by the coordinate transformation. 558 The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}. 559 \item[\np{nn\_bathy}~\forcode{= -1}]: 560 a domain with a bump of topography one third of the domain width at the central latitude. 561 This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount. 562 \item[\np{nn\_bathy}~\forcode{= 1}]: 563 read a bathymetry and ice shelf draft (if needed). 564 The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at 565 each grid point of the model grid. 566 The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto 567 the horizontal ocean mesh. 568 Defining the bathymetry also defines the coastline: where the bathymetry is zero, 569 no model levels are defined (all levels are masked). 570 571 The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at 572 each grid point of the model grid. 573 This file is only needed if \np{ln\_isfcav}~\forcode{= .true.}. 574 Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. 575 \end{description} 576 577 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies 578 (\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with 579 the rest of the ocean. 610 611 Note that, without ice shelves cavities, 612 masks at $t-$ and $w-$points are identical with the numerical indexing used 613 (\autoref{subsec:DOM_Num_Index}). 614 Nevertheless, 615 $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 616 exactly in the same way as for the bottom boundary. 617 618 %% The specification of closed lateral boundaries requires that at least 619 %% the first and last rows and columns of the \textit{mbathy} array are set to zero. 620 %% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 621 %% the second one and its first column equal to the last but one (and so too the mask arrays) 622 %% (see \autoref{fig:LBC_jperio}). 623 624 % Closed seas 625 %% ================================================================================================= 626 \subsection{Closed seas} 627 \label{subsec:DOM_closea} 628 629 When a global ocean is coupled to an atmospheric model it is better to 630 represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if 631 the model resolution does not allow their communication with the rest of the ocean. 580 632 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 581 633 so these seas can be removed from the ocean domain. 582 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}), 583 but the code has to be adapted to the user's configuration. 584 585 % ------------------------------------------------------------------------------------------------------------- 586 % z-coordinate and reference coordinate transformation 587 % ------------------------------------------------------------------------------------------------------------- 588 \subsection[$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and ref. coordinate] 589 {$Z$-coordinate (\protect\np{ln\_zco}~\forcode{= .true.}) and reference coordinate} 590 \label{subsec:DOM_zco} 591 592 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 593 \begin{figure}[!tb] 594 \begin{center} 595 \includegraphics[]{Fig_zgr} 596 \caption{ 597 \protect\label{fig:zgr} 598 Default vertical mesh for ORCA2: 30 ocean levels (L30). 599 Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from 600 \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate. 601 } 602 \end{center} 603 \end{figure} 604 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 605 606 The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points, 607 respectively. 608 As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels. 609 $gdepw_0(1)$ is the ocean surface. 610 There are at most \jp{jpk}-1 $t$-points inside the ocean, 611 the additional $t$-point at $jk = jpk$ is below the sea floor and is not used. 612 The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose 613 analytical derivative with respect to $k$ provides the vertical scale factors. 614 The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$. 615 This is done in routine \mdl{domzgr} through statement functions, 616 using parameters provided in the \ngn{namcfg} namelist. 617 618 It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}~\forcode{= 0}). 619 In that case, the parameters \jp{jpk} (number of $w$-levels) and 620 \np{pphmax} (total ocean depth in meters) fully define the grid. 621 622 For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface. 623 The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps): 624 \begin{gather} 625 \label{eq:DOM_zgr_ana_1} 626 z_0 (k) = h_{sur} - h_0 \; k - \; h_1 \; \log \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\ 627 e_3^0(k) = \lt| - h_0 - h_1 \; \tanh \big[ (k - h_{th}) / h_{cr} \big] \rt| 628 \end{gather} 629 where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels. 630 Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with 631 a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}). 632 633 If the ice shelf cavities are opened (\np{ln\_isfcav}~\forcode{= .true.}), the definition of $z_0$ is the same. 634 However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: 635 \begin{equation} 636 \label{eq:DOM_zgr_ana_2} 637 \begin{split} 638 e_3^T(k) &= z_W (k + 1) - z_W (k ) \\ 639 e_3^W(k) &= z_T (k ) - z_T (k - 1) 640 \end{split} 641 \end{equation} 642 This formulation decrease the self-generated circulation into the ice shelf cavity 643 (which can, in extreme case, leads to blow up).\\ 644 645 The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and 646 a depth which varies from 0 at the sea surface to a minimum of $-5000~m$. 647 This leads to the following conditions: 648 \begin{equation} 649 \label{eq:DOM_zgr_coef} 650 \begin{array}{ll} 651 e_3 (1 + 1/2) = 10. & z(1 ) = 0. \\ 652 e_3 (jpk - 1/2) = 500. & z(jpk) = -5000. 653 \end{array} 654 \end{equation} 655 656 With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$, 657 the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in 658 \autoref{eq:DOM_zgr_ana_2} have been determined such that 659 \autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method. 660 For the first standard ORCA2 vertical grid this led to the following values: 661 $h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$. 662 The resulting depths and scale factors as a function of the model levels are shown in 663 \autoref{fig:zgr} and given in \autoref{tab:orca_zgr}. 664 Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. 665 666 Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them. 667 In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$., 668 in \ngn{namcfg} namelist, and specifies instead the four following parameters: 669 \begin{itemize} 670 \item 671 \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional). 672 The larger \np{ppacr}, the smaller the stretching. 673 Values from $3$ to $10$ are usual. 674 \item 675 \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs 676 (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) 677 \item 678 \np{ppdzmin}: minimum thickness for the top layer (in meters). 679 \item 680 \np{pphmax}: total depth of the ocean (meters). 681 \end{itemize} 682 As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are: 683 \jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$. 684 685 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 686 \begin{table} 687 \begin{center} 688 \begin{tabular}{c||r|r|r|r} 689 \hline 690 \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\ 691 \hline 692 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ 693 \hline 694 2 & \textbf{ 15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ 695 \hline 696 3 & \textbf{ 25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ 697 \hline 698 4 & \textbf{ 35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ 699 \hline 700 5 & \textbf{ 45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ 701 \hline 702 6 & \textbf{ 55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ 703 \hline 704 7 & \textbf{ 65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ 705 \hline 706 8 & \textbf{ 75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ 707 \hline 708 9 & \textbf{ 85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ 709 \hline 710 10 & \textbf{ 95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ 711 \hline 712 11 & \textbf{ 105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ 713 \hline 714 12 & \textbf{ 116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ 715 \hline 716 13 & \textbf{ 128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ 717 \hline 718 14 & \textbf{ 142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ 719 \hline 720 15 & \textbf{ 158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ 721 \hline 722 16 & \textbf{ 181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ 723 \hline 724 17 & \textbf{ 216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ 725 \hline 726 18 & \textbf{ 272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ 727 \hline 728 19 & \textbf{ 364.30} & 312.74 & \textbf{ 116.11} & 90.99 \\ 729 \hline 730 20 & \textbf{ 511.53} & 429.72 & \textbf{ 181.55} & 146.43 \\ 731 \hline 732 21 & \textbf{ 732.20} & 611.89 & \textbf{ 261.03} & 220.35 \\ 733 \hline 734 22 & \textbf{ 1033.22} & 872.87 & \textbf{ 339.39} & 301.42 \\ 735 \hline 736 23 & \textbf{ 1405.70} & 1211.59 & \textbf{ 402.26} & 373.31 \\ 737 \hline 738 24 & \textbf{ 1830.89} & 1612.98 & \textbf{ 444.87} & 426.00 \\ 739 \hline 740 25 & \textbf{ 2289.77} & 2057.13 & \textbf{ 470.55} & 459.47 \\ 741 \hline 742 26 & \textbf{ 2768.24} & 2527.22 & \textbf{ 484.95} & 478.83 \\ 743 \hline 744 27 & \textbf{ 3257.48} & 3011.90 & \textbf{ 492.70} & 489.44 \\ 745 \hline 746 28 & \textbf{ 3752.44} & 3504.46 & \textbf{ 496.78} & 495.07 \\ 747 \hline 748 29 & \textbf{ 4250.40} & 4001.16 & \textbf{ 498.90} & 498.02 \\ 749 \hline 750 30 & \textbf{ 4749.91} & 4500.02 & \textbf{ 500.00} & 499.54 \\ 751 \hline 752 31 & \textbf{ 5250.23} & 5000.00 & \textbf{ 500.56} & 500.33 \\ 753 \hline 754 \end{tabular} 755 \end{center} 756 \caption{ 757 \protect\label{tab:orca_zgr} 758 Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from 759 \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef} 760 } 761 \end{table} 762 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 763 764 % ------------------------------------------------------------------------------------------------------------- 765 % z-coordinate with partial step 766 % ------------------------------------------------------------------------------------------------------------- 767 \subsection{$Z$-coordinate with partial step (\protect\np{ln\_zps}~\forcode{= .true.})} 768 \label{subsec:DOM_zps} 769 %--------------------------------------------namdom------------------------------------------------------- 770 771 \nlst{namdom} 772 %-------------------------------------------------------------------------------------------------------------- 773 774 In $z$-coordinate partial step, 775 the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in 776 the previous section, \textit{except} in the bottom layer. 777 The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to 778 allow a better representation of the bathymetry, especially in the case of small slopes 779 (where the bathymetry varies by less than one level thickness from one grid point to the next). 780 The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry. 781 With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. 782 The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: 783 the maximum thickness allowed is $2*e_{3t}(jpk - 1)$. 784 This has to be kept in mind when specifying values in \ngn{namdom} namelist, 785 as the maximum depth \np{pphmax} in partial steps: 786 for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid, 787 the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$). 788 Two variables in the namdom namelist are used to define the partial step vertical grid. 789 The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is 790 the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat} 791 (a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$). 792 793 \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } } 794 795 % ------------------------------------------------------------------------------------------------------------- 796 % s-coordinate 797 % ------------------------------------------------------------------------------------------------------------- 798 \subsection{$S$-coordinate (\protect\np{ln\_sco}~\forcode{= .true.})} 799 \label{subsec:DOM_sco} 800 %------------------------------------------nam_zgr_sco--------------------------------------------------- 801 % 802 %\nlst{namzgr_sco} 803 %-------------------------------------------------------------------------------------------------------------- 804 Options are defined in \ngn{namzgr\_sco}. 805 In $s$-coordinate (\np{ln\_sco}~\forcode{= .true.}), the depth and thickness of the model levels are defined from 806 the product of a depth field and either a stretching function or its derivative, respectively: 807 808 \begin{align*} 809 % \label{eq:DOM_sco_ana} 810 z(k) &= h(i,j) \; z_0 (k) \\ 811 e_3(k) &= h(i,j) \; z_0'(k) 812 \end{align*} 813 814 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and 815 $z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom. 816 The depth field $h$ is not necessary the ocean depth, 817 since a mixed step-like and bottom-following representation of the topography can be used 818 (\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}). 819 The namelist parameter \np{rn\_rmax} determines the slope at which 820 the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. 821 The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as 822 the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 823 824 Options for stretching the coordinate are provided as examples, 825 but care must be taken to ensure that the vertical stretch used is appropriate for the application. 826 827 The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true 828 (\np{ln\_s\_SH94}~\forcode{= .false.} and \np{ln\_s\_SF12}~\forcode{= .false.}). 829 This uses a depth independent $\tanh$ function for the stretching \citep{madec.delecluse.ea_JPO96}: 830 831 \[ 832 z = s_{min} + C (s) (H - s_{min}) 833 % \label{eq:SH94_1} 834 \] 835 836 where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 837 allows a $z$-coordinate to placed on top of the stretched coordinate, 838 and $z$ is the depth (negative down from the asea surface). 839 \begin{gather*} 840 s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1 841 % \label{eq:DOM_s} 842 \\ 843 % \label{eq:DOM_sco_function} 844 C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)} 845 \end{gather*} 846 847 A stretching function, 848 modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np{ln\_s\_SH94}~\forcode{= .true.}), 849 is also available and is more commonly used for shelf seas modelling: 850 851 \[ 852 C(s) = (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)} 853 + b \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] - \tanh \lt( \frac{\theta}{2} \rt)} 854 { 2 \tanh \lt( \frac{\theta}{2} \rt)} 855 % \label{eq:SH94_2} 856 \] 857 858 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 859 \begin{figure}[!ht] 860 \begin{center} 861 \includegraphics[]{Fig_sco_function} 862 \caption{ 863 \protect\label{fig:sco_function} 864 Examples of the stretching function applied to a seamount; 865 from left to right: surface, surface and bottom, and bottom intensified resolutions 866 } 867 \end{center} 868 \end{figure} 869 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 870 871 where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to 872 the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 873 bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$. 874 $b$ has been designed to allow surface and/or bottom increase of the vertical resolution 875 (\autoref{fig:sco_function}). 876 877 Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in 878 an analytical terrain-following stretching \citet{siddorn.furner_OM13}. 879 In this case the a stretching function $\gamma$ is defined such that: 880 881 \begin{equation} 882 z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1 883 % \label{eq:z} 884 \end{equation} 885 886 The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: 887 888 \begin{gather*} 889 % \label{eq:DOM_gamma_deriv} 890 \gamma = A \lt( \sigma - \frac{1}{2} (\sigma^2 + f (\sigma)) \rt) 891 + B \lt( \sigma^3 - f (\sigma) \rt) + f (\sigma) \\ 892 \intertext{Where:} 893 % \label{eq:DOM_gamma} 894 f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2} 895 \quad \text{and} \quad \sigma = \frac{k}{n - 1} 896 \end{gather*} 897 898 This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of 899 the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards 900 the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and 901 user prescribed surface (\np{rn\_zs}) and bottom depths. 902 The bottom cell depth in this example is given as a function of water depth: 903 904 \[ 905 % \label{eq:DOM_zb} 906 Z_b = h a + b 907 \] 908 909 where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. 910 911 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 912 \begin{figure}[!ht] 913 \includegraphics[]{Fig_DOM_compare_coordinates_surface} 914 \caption{ 915 A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines), 916 a 50 level $Z$-coordinate (contoured surfaces) and 917 the \citet{siddorn.furner_OM13} $S$-coordinate (dashed lines) in the surface $100~m$ for 918 a idealised bathymetry that goes from $50~m$ to $5500~m$ depth. 919 For clarity every third coordinate surface is shown. 920 } 921 \label{fig:fig_compare_coordinates_surface} 922 \end{figure} 923 % >>>>>>>>>>>>>>>>>>>>>>>>>>>> 924 925 This gives a smooth analytical stretching in computational space that is constrained to 926 given specified surface and bottom grid cell thicknesses in real space. 927 This is not to be confused with the hybrid schemes that 928 superimpose geopotential coordinates on terrain following coordinates thus 929 creating a non-analytical vertical coordinate that 930 therefore may suffer from large gradients in the vertical resolutions. 931 This stretching is less straightforward to implement than the \citet{song.haidvogel_JCP94} stretching, 932 but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. 933 934 As with the \citet{song.haidvogel_JCP94} stretching the stretch is only applied at depths greater than 935 the critical depth $h_c$. 936 In this example two options are available in depths shallower than $h_c$, 937 with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false 938 (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$). 939 940 Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as 941 large slopes lead to hydrostatic consistency. 942 A hydrostatic consistency parameter diagnostic following \citet{haney_JPO91} has been implemented, 943 and is output as part of the model mesh file at the start of the run. 944 945 % ------------------------------------------------------------------------------------------------------------- 946 % z*- or s*-coordinate 947 % ------------------------------------------------------------------------------------------------------------- 948 \subsection{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}~\forcode{= .false.})} 949 \label{subsec:DOM_zgr_star} 950 951 This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site. 952 953 %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances 954 955 % ------------------------------------------------------------------------------------------------------------- 956 % level bathymetry and mask 957 % ------------------------------------------------------------------------------------------------------------- 958 \subsection{Level bathymetry and mask} 959 \label{subsec:DOM_msk} 960 961 Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with 962 steps that follow the face of the model cells (step like topography) \citep{madec.delecluse.ea_JPO96}. 963 The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which 964 gives the number of ocean levels (\ie those that are not masked) at each $t$-point. 965 mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which 966 gdept $\leq$ bathy. 967 968 Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after 969 mbathy is computed. 970 Isolated grid points that do not communicate with another ocean point at the same level are eliminated. 971 972 As for the representation of bathymetry, a 2D integer array, misfdep, is created. 973 misfdep defines the level of the first wet $t$-point. 974 All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked. 975 By default, $misfdep(:,:) = 1$ and no cells are masked. 976 977 In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 978 the cavities are performed in the \textit{zgr\_isf} routine. 979 The compatibility between ice shelf draft and bathymetry is checked. 980 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked). 981 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, 982 the bathymetry or the ice shelf draft is dug to fit this constrain. 983 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked. 984 985 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 986 \begin{alignat*}{2} 987 tmask(i,j,k) &= & & 988 \begin{cases} 989 0 &\text{if $ k < misfdep(i,j)$} \\ 990 1 &\text{if $misfdep(i,j) \leq k \leq mbathy(i,j)$} \\ 991 0 &\text{if $ k > mbathy(i,j)$} 992 \end{cases} 993 \\ 994 umask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 995 vmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j + 1,k) \\ 996 fmask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 997 & &* &tmask(i,j,k) * tmask(i + 1,j, k) \\ 998 wmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j,k - 1) \\ 999 \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 1000 \end{alignat*} 1001 1002 Note that, without ice shelves cavities, 1003 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 1004 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 1005 exactly in the same way as for the bottom boundary. 1006 1007 The specification of closed lateral boundaries requires that at least 1008 the first and last rows and columns of the \textit{mbathy} array are set to zero. 1009 In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to 1010 the second one and its first column equal to the last but one (and so too the mask arrays) 1011 (see \autoref{fig:LBC_jperio}). 1012 1013 % ================================================================ 1014 % Domain: Initial State (dtatsd & istate) 1015 % ================================================================ 1016 \section{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 1017 \label{sec:DTA_tsd} 1018 %-----------------------------------------namtsd------------------------------------------- 1019 1020 \nlst{namtsd} 1021 %------------------------------------------------------------------------------------------ 1022 1023 Options are defined in \ngn{namtsd}. 1024 By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and 1025 salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. 634 The user has the option to 635 set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to 636 optionally decide on the fate of any freshwater imbalance over the area. 637 The options are explained in \autoref{sec:MISC_closea} but 638 it should be noted here that a successful use of these options requires 639 appropriate mask fields to be present in the domain configuration file. 640 Among the possibilities are: 641 642 \begin{clines} 643 int closea_mask /* non-zero values in closed sea areas for optional masking */ 644 int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ 645 int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp) */ 646 \end{clines} 647 648 %% ================================================================================================= 649 \subsection{Output grid files} 650 \label{subsec:DOM_meshmask} 651 652 Most of the arrays relating to a particular ocean model configuration discussed in this chapter 653 (grid-point position, scale factors) can be saved in a file if 654 namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to 655 \forcode{.true.}; 656 the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. 657 This is only really useful if 658 the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and 659 checking or confirmation is required. 660 661 Alternatively, all the arrays relating to a particular ocean model configuration 662 (grid-point position, scale factors, depths and masks) can be saved in 663 a file called \texttt{mesh\_mask} if 664 namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to 665 \forcode{.true.}. 666 This file contains additional fields that can be useful for post-processing applications. 667 668 %% ================================================================================================= 669 \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 670 \label{sec:DOM_DTA_tsd} 671 672 \begin{listing} 673 \nlst{namtsd} 674 \caption{\forcode{&namtsd}} 675 \label{lst:namtsd} 676 \end{listing} 677 678 Basic initial state options are defined in \nam{tsd}{tsd}. 679 By default, the ocean starts from rest (the velocity field is set to zero) and 680 the initialization of temperature and salinity fields is controlled through the \np{ln_tsd_init}{ln\_tsd\_init} namelist parameter. 681 1026 682 \begin{description} 1027 \item [\np{ln\_tsd\_init}~\forcode{= .true.}]1028 use a T and S input files that can be given on the model grid itself or on their native input data grid.683 \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on 684 the model grid itself or on their native input data grids. 1029 685 In the latter case, 1030 686 the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 1031 687 (see \autoref{subsec:SBC_iof}). 1032 The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. 688 The information relating to the input files are specified in 689 the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 1033 690 The computation is done in the \mdl{dtatsd} module. 1034 \item[\np{ln\_tsd\_init}~\forcode{= .false.}] 1035 use constant salinity value of $35.5~psu$ and an analytical profile of temperature 1036 (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. 691 \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via 692 a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 693 The default version sets horizontally uniform T and profiles as used in the GYRE configuration 694 (see \autoref{sec:CFGS_gyre}). 1037 695 \end{description} 1038 696 1039 \biblio 1040 1041 \pindex 697 \onlyinsubfile{\input{../../global/epilogue}} 1042 698 1043 699 \end{document}
Note: See TracChangeset
for help on using the changeset viewer.