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