1 | % ================================================================ |
---|
2 | % Chapter 2 ——— Space and Time Domain (DOM) |
---|
3 | % ================================================================ |
---|
4 | \chapter{Space Domain (DOM) } |
---|
5 | \label{DOM} |
---|
6 | \minitoc |
---|
7 | |
---|
8 | % Missing things: |
---|
9 | % - istate: description of the initial state ==> this has to be put elsewhere.. |
---|
10 | % perhaps in MISC ? By the way the initialisation of T S and dynamics |
---|
11 | % should be put outside of DOM routine (better with TRC staff and off-line |
---|
12 | % tracers) |
---|
13 | % -geo2ocean: how to switch from geographic to mesh coordinate |
---|
14 | % - domclo: closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled |
---|
15 | |
---|
16 | |
---|
17 | \newpage |
---|
18 | $\ $\newline % force a new ligne |
---|
19 | |
---|
20 | Having defined the continuous equations in Chap.~\ref{PE} and chosen a time |
---|
21 | discretization Chap.~\ref{STP}, we need to choose a discretization on a grid, |
---|
22 | and numerical algorithms. In the present chapter, we provide a general description |
---|
23 | of the staggered grid used in \NEMO, and other information relevant to the main |
---|
24 | directory routines as well as the DOM (DOMain) directory. |
---|
25 | |
---|
26 | $\ $\newline % force a new lign |
---|
27 | |
---|
28 | % ================================================================ |
---|
29 | % Fundamentals of the Discretisation |
---|
30 | % ================================================================ |
---|
31 | \section{Fundamentals of the Discretisation} |
---|
32 | \label{DOM_basics} |
---|
33 | |
---|
34 | % ------------------------------------------------------------------------------------------------------------- |
---|
35 | % Arrangement of Variables |
---|
36 | % ------------------------------------------------------------------------------------------------------------- |
---|
37 | \subsection{Arrangement of Variables} |
---|
38 | \label{DOM_cell} |
---|
39 | |
---|
40 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
41 | \begin{figure}[!tb] \begin{center} |
---|
42 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} |
---|
43 | \caption{ \label{Fig_cell} |
---|
44 | Arrangement of variables. $t$ indicates scalar points where temperature, |
---|
45 | salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) |
---|
46 | indicates vector points, and $f$ indicates vorticity points where both relative and |
---|
47 | planetary vorticities are defined} |
---|
48 | \end{center} \end{figure} |
---|
49 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
50 | |
---|
51 | The numerical techniques used to solve the Primitive Equations in this model are |
---|
52 | based on the traditional, centred second-order finite difference approximation. |
---|
53 | Special attention has been given to the homogeneity of the solution in the three |
---|
54 | space directions. The arrangement of variables is the same in all directions. |
---|
55 | It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector |
---|
56 | points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). |
---|
57 | This is the generalisation to three dimensions of the well-known ``C'' grid in |
---|
58 | Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and |
---|
59 | planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge |
---|
60 | and the barotropic stream function $\psi$ is defined at horizontal points overlying |
---|
61 | the $\zeta$ and $f$-points. |
---|
62 | |
---|
63 | The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined |
---|
64 | by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. |
---|
65 | The grid-points are located at integer or integer and a half value of $(i,j,k)$ as |
---|
66 | indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$, |
---|
67 | $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale |
---|
68 | factors are defined. Each scale factor is defined as the local analytical value |
---|
69 | provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial |
---|
70 | derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and |
---|
71 | $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. |
---|
72 | Discrete partial derivatives are formulated by the traditional, centred second order |
---|
73 | finite difference approximation while the scale factors are chosen equal to their |
---|
74 | local analytical value. An important point here is that the partial derivative of the |
---|
75 | scale factors must be evaluated by centred finite difference approximation, not |
---|
76 | from their analytical expression. This preserves the symmetry of the discrete set |
---|
77 | of equations and therefore satisfies many of the continuous properties (see |
---|
78 | Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain |
---|
79 | size: when needed, an area, volume, or the total ocean depth must be evaluated |
---|
80 | as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). |
---|
81 | |
---|
82 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
83 | \begin{table}[!tb] |
---|
84 | \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} |
---|
85 | \hline |
---|
86 | T &$i$ & $j$ & $k$ \\ \hline |
---|
87 | u & $i+1/2$ & $j$ & $k$ \\ \hline |
---|
88 | v & $i$ & $j+1/2$ & $k$ \\ \hline |
---|
89 | w & $i$ & $j$ & $k+1/2$ \\ \hline |
---|
90 | f & $i+1/2$ & $j+1/2$ & $k$ \\ \hline |
---|
91 | uw & $i+1/2$ & $j$ & $k+1/2$ \\ \hline |
---|
92 | vw & $i$ & $j+1/2$ & $k+1/2$ \\ \hline |
---|
93 | fw & $i+1/2$ & $j+1/2$ & $k+1/2$ \\ \hline |
---|
94 | \end{tabular} |
---|
95 | \caption{ \label{Tab_cell} |
---|
96 | Location of grid-points as a function of integer or integer and a half value of the column, |
---|
97 | line or level. This indexing is only used for the writing of the semi-discrete equation. |
---|
98 | In the code, the indexing uses integer values only and has a reverse direction |
---|
99 | in the vertical (see \S\ref{DOM_Num_Index})} |
---|
100 | \end{center} |
---|
101 | \end{table} |
---|
102 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
103 | |
---|
104 | % ------------------------------------------------------------------------------------------------------------- |
---|
105 | % Vector Invariant Formulation |
---|
106 | % ------------------------------------------------------------------------------------------------------------- |
---|
107 | \subsection{Discrete Operators} |
---|
108 | \label{DOM_operators} |
---|
109 | |
---|
110 | Given the values of a variable $q$ at adjacent points, the differencing and |
---|
111 | averaging operators at the midpoint between them are: |
---|
112 | \begin{subequations} \label{Eq_di_mi} |
---|
113 | \begin{align} |
---|
114 | \delta _i [q] &= \ \ q(i+1/2) - q(i-1/2) \\ |
---|
115 | \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2 |
---|
116 | \end{align} |
---|
117 | \end{subequations} |
---|
118 | |
---|
119 | Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and |
---|
120 | $k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a |
---|
121 | variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- |
---|
122 | and $w$-points while its Laplacien is defined at $t$-point. These operators have |
---|
123 | the following discrete forms in the curvilinear $s$-coordinate system: |
---|
124 | \begin{equation} \label{Eq_DOM_grad} |
---|
125 | \nabla q\equiv \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,\mathbf{i} |
---|
126 | + \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,\mathbf{j} |
---|
127 | + \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,\mathbf{k} |
---|
128 | \end{equation} |
---|
129 | \begin{multline} \label{Eq_DOM_lap} |
---|
130 | \Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} } |
---|
131 | \;\left( \delta_i \left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right] |
---|
132 | + \delta_j \left[ \frac{e_{1v}\,e_{3v}} {e_{2v}} \;\delta_{j+1/2} [q] \right] \; \right) \\ |
---|
133 | +\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} } \;\delta_{k+1/2} [q] \right] |
---|
134 | \end{multline} |
---|
135 | |
---|
136 | Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ |
---|
137 | defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, |
---|
138 | and $f$-points, and its divergence defined at $t$-points: |
---|
139 | \begin{eqnarray} \label{Eq_DOM_curl} |
---|
140 | \nabla \times {\rm{\bf A}}\equiv & |
---|
141 | \frac{1}{e_{2v}\,e_{3vw} } \ \left( \delta_{j +1/2} \left[e_{3w}\,a_3 \right] -\delta_{k+1/2} \left[e_{2v} \,a_2 \right] \right) &\ \mathbf{i} \\ |
---|
142 | +& \frac{1}{e_{2u}\,e_{3uw}} \ \left( \delta_{k+1/2} \left[e_{1u}\,a_1 \right] -\delta_{i +1/2} \left[e_{3w}\,a_3 \right] \right) &\ \mathbf{j} \\ |
---|
143 | +& \frac{1}{e_{1f} \,e_{2f} } \ \left( \delta_{i +1/2} \left[e_{2v}\,a_2 \right] -\delta_{j +1/2} \left[e_{1u}\,a_1 \right] \right) &\ \mathbf{k} |
---|
144 | \end{eqnarray} |
---|
145 | \begin{equation} \label{Eq_DOM_div} |
---|
146 | \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right] |
---|
147 | +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right] |
---|
148 | \end{equation} |
---|
149 | |
---|
150 | In the special case of a pure $z$-coordinate system, \eqref{Eq_DOM_lap} and |
---|
151 | \eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor |
---|
152 | becomes a function of the single variable $k$ and thus does not depend on the |
---|
153 | horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to: |
---|
154 | \begin{equation*} |
---|
155 | \nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right] |
---|
156 | +\delta_j \left[e_{1v}\, a_2 \right] \right) |
---|
157 | +\frac{1}{e_{3t}} \delta_k \left[ a_3 \right] |
---|
158 | \end{equation*} |
---|
159 | |
---|
160 | The vertical average over the whole water column denoted by an overbar becomes |
---|
161 | for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): |
---|
162 | \begin{equation} \label{DOM_bar} |
---|
163 | \bar q = \frac{1}{H} \int_{k^b}^{k^o} {q\;e_{3q} \,dk} |
---|
164 | \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } |
---|
165 | \end{equation} |
---|
166 | where $H_q$ is the ocean depth, which is the masked sum of the vertical scale |
---|
167 | factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, |
---|
168 | and the symbol $k^o$ refers to a summation over all grid points of the same type |
---|
169 | in the direction indicated by the subscript (here $k$). |
---|
170 | |
---|
171 | In continuous form, the following properties are satisfied: |
---|
172 | \begin{equation} \label{Eq_DOM_curl_grad} |
---|
173 | \nabla \times \nabla q ={\rm {\bf {0}}} |
---|
174 | \end{equation} |
---|
175 | \begin{equation} \label{Eq_DOM_div_curl} |
---|
176 | \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0 |
---|
177 | \end{equation} |
---|
178 | |
---|
179 | It is straightforward to demonstrate that these properties are verified locally in |
---|
180 | discrete form as soon as the scalar $q$ is taken at $t$-points and the vector |
---|
181 | \textbf{A} has its components defined at vector points $(u,v,w)$. |
---|
182 | |
---|
183 | Let $a$ and $b$ be two fields defined on the mesh, with value zero inside |
---|
184 | continental area. Using integration by parts it can be shown that the differencing |
---|
185 | operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, |
---|
186 | and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, |
---|
187 | $\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear |
---|
188 | operators, $i.e.$ |
---|
189 | \begin{align} |
---|
190 | \label{DOM_di_adj} |
---|
191 | \sum\limits_i { a_i \;\delta _i \left[ b \right]} |
---|
192 | &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} } \\ |
---|
193 | \label{DOM_mi_adj} |
---|
194 | \sum\limits_i { a_i \;\overline b^{\,i}} |
---|
195 | & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } |
---|
196 | \end{align} |
---|
197 | |
---|
198 | In other words, the adjoint of the differencing and averaging operators are |
---|
199 | $\delta_i^*=\delta_{i+1/2}$ and |
---|
200 | ${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. |
---|
201 | These two properties will be used extensively in the Appendix~\ref{Apdx_C} to |
---|
202 | demonstrate integral conservative properties of the discrete formulation chosen. |
---|
203 | |
---|
204 | % ------------------------------------------------------------------------------------------------------------- |
---|
205 | % Numerical Indexing |
---|
206 | % ------------------------------------------------------------------------------------------------------------- |
---|
207 | \subsection{Numerical Indexing} |
---|
208 | \label{DOM_Num_Index} |
---|
209 | |
---|
210 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
211 | \begin{figure}[!tb] \begin{center} |
---|
212 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf} |
---|
213 | \caption{ \label{Fig_index_hor} |
---|
214 | Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates |
---|
215 | the cell in which variables contained in arrays have the same $i$- and $j$-indices} |
---|
216 | \end{center} \end{figure} |
---|
217 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
218 | |
---|
219 | The array representation used in the \textsc{Fortran} code requires an integer |
---|
220 | indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is |
---|
221 | associated with the use of integer values for $t$-points and both integer and |
---|
222 | integer and a half values for all the other points. Therefore a specific integer |
---|
223 | indexing must be defined for points other than $t$-points ($i.e.$ velocity and |
---|
224 | vorticity grid-points). Furthermore, the direction of the vertical indexing has |
---|
225 | been changed so that the surface level is at $k=1$. |
---|
226 | |
---|
227 | % ----------------------------------- |
---|
228 | % Horizontal Indexing |
---|
229 | % ----------------------------------- |
---|
230 | \subsubsection{Horizontal Indexing} |
---|
231 | \label{DOM_Num_Index_hor} |
---|
232 | |
---|
233 | The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. |
---|
234 | For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point |
---|
235 | (northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}). |
---|
236 | A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. |
---|
237 | |
---|
238 | % ----------------------------------- |
---|
239 | % Vertical indexing |
---|
240 | % ----------------------------------- |
---|
241 | \subsubsection{Vertical Indexing} |
---|
242 | \label{DOM_Num_Index_vertical} |
---|
243 | |
---|
244 | In the vertical, the chosen indexing requires special attention since the |
---|
245 | $k$-axis is re-orientated downward in the \textsc{Fortran} code compared |
---|
246 | to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}. |
---|
247 | The sea surface corresponds to the $w$-level $k=1$ which is the same index |
---|
248 | as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) |
---|
249 | either corresponds to the ocean floor or is inside the bathymetry while the last |
---|
250 | $t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that |
---|
251 | for an increasing $k$ index, a $w$-point and the $t$-point just below have the |
---|
252 | same $k$ index, in opposition to what is done in the horizontal plane where |
---|
253 | it is the $t$-point and the nearest velocity points in the direction of the horizontal |
---|
254 | axis that have the same $i$ or $j$ index (compare the dashed area in |
---|
255 | Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are |
---|
256 | chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} |
---|
257 | code \emph{before all the vertical derivatives} of the discrete equations given in |
---|
258 | this documentation. |
---|
259 | |
---|
260 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
261 | \begin{figure}[!pt] \begin{center} |
---|
262 | \includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf} |
---|
263 | \caption{ \label{Fig_index_vert} |
---|
264 | Vertical integer indexing used in the \textsc{Fortran } code. Note that |
---|
265 | the $k$-axis is orientated downward. The dashed area indicates the cell in |
---|
266 | which variables contained in arrays have the same $k$-index.} |
---|
267 | \end{center} \end{figure} |
---|
268 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
269 | |
---|
270 | % ----------------------------------- |
---|
271 | % Domain Size |
---|
272 | % ----------------------------------- |
---|
273 | \subsubsection{Domain Size} |
---|
274 | \label{DOM_size} |
---|
275 | |
---|
276 | The total size of the computational domain is set by the parameters \np{jpiglo}, |
---|
277 | \np{jpjglo} and \np{jpkdta} in the $i$, $j$ and $k$ directions respectively. They are |
---|
278 | given as namelist variables in the \ngn{namcfg} namelist. |
---|
279 | |
---|
280 | Note that are other namelist variables in the \ngn{namcfg} namelist that refer to |
---|
281 | the domain size. |
---|
282 | The two variables \np{jpidta} and \np{jpjdta} may be larger than \np{jpiglo}, \np{jpjglo} |
---|
283 | when the user wants to use only a sub-region of a given configuration. This is |
---|
284 | the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of |
---|
285 | the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters |
---|
286 | $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is |
---|
287 | run in parallel using domain decomposition (\key{mpp\_mpi} defined, see |
---|
288 | \S\ref{LBC_mpp}). |
---|
289 | |
---|
290 | |
---|
291 | $\ $\newline % force a new lign |
---|
292 | |
---|
293 | % ================================================================ |
---|
294 | % Domain: Horizontal Grid (mesh) |
---|
295 | % ================================================================ |
---|
296 | \section [Domain: Horizontal Grid (mesh) (\textit{domhgr})] |
---|
297 | {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} } |
---|
298 | \label{DOM_hgr} |
---|
299 | |
---|
300 | % ------------------------------------------------------------------------------------------------------------- |
---|
301 | % Coordinates and scale factors |
---|
302 | % ------------------------------------------------------------------------------------------------------------- |
---|
303 | \subsection{Coordinates and scale factors} |
---|
304 | \label{DOM_hgr_coord_e} |
---|
305 | |
---|
306 | The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined |
---|
307 | by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. |
---|
308 | The grid-points are located at integer or integer and a half values of as indicated |
---|
309 | in Table~\ref{Tab_cell}. The associated scale factors are defined using the |
---|
310 | analytical first derivative of the transformation \eqref{Eq_scale_factors}. These |
---|
311 | definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which |
---|
312 | provide the horizontal and vertical meshes, respectively. This section deals with |
---|
313 | the horizontal mesh parameters. |
---|
314 | |
---|
315 | In a horizontal plane, the location of all the model grid points is defined from the |
---|
316 | analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a |
---|
317 | function of $(i,j)$. The horizontal scale factors are calculated using |
---|
318 | \eqref{Eq_scale_factors}. For example, when the longitude and latitude are |
---|
319 | function of a single value ($i$ and $j$, respectively) (geographical configuration |
---|
320 | of the mesh), the horizontal mesh definition reduces to define the wanted |
---|
321 | $\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the |
---|
322 | \mdl{domhgr} module. The model computes the grid-point positions and scale |
---|
323 | factors in the horizontal plane as follows: |
---|
324 | \begin{flalign*} |
---|
325 | \lambda_t &\equiv \text{glamt}= \lambda(i) & \varphi_t &\equiv \text{gphit} = \varphi(j)\\ |
---|
326 | \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ |
---|
327 | \lambda_v &\equiv \text{glamv}= \lambda(i) & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ |
---|
328 | \lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) |
---|
329 | \end{flalign*} |
---|
330 | \begin{flalign*} |
---|
331 | e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j) |& |
---|
332 | e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)| \\ |
---|
333 | e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2) \; \cos\varphi(j) |& |
---|
334 | e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ |
---|
335 | e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j+1/2) |& |
---|
336 | e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\ |
---|
337 | e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2) |& |
---|
338 | e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)| |
---|
339 | \end{flalign*} |
---|
340 | where the last letter of each computational name indicates the grid point |
---|
341 | considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with |
---|
342 | all universal constants). Note that the horizontal position of and scale factors |
---|
343 | at $w$-points are exactly equal to those of $t$-points, thus no specific arrays |
---|
344 | are defined at $w$-points. |
---|
345 | |
---|
346 | Note that the definition of the scale factors ($i.e.$ as the analytical first derivative |
---|
347 | of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is |
---|
348 | specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined |
---|
349 | locally at a $t$-point, whereas many other models on a C grid choose to define |
---|
350 | such a scale factor as the distance between the $U$-points on each side of the |
---|
351 | $t$-point. Relying on an analytical transformation has two advantages: firstly, there |
---|
352 | is no ambiguity in the scale factors appearing in the discrete equations, since they |
---|
353 | are first introduced in the continuous equations; secondly, analytical transformations |
---|
354 | encourage good practice by the definition of smoothly varying grids (rather than |
---|
355 | allowing the user to set arbitrary jumps in thickness between adjacent layers) |
---|
356 | \citep{Treguier1996}. An example of the effect of such a choice is shown in |
---|
357 | Fig.~\ref{Fig_zgr_e3}. |
---|
358 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
359 | \begin{figure}[!t] \begin{center} |
---|
360 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf} |
---|
361 | \caption{ \label{Fig_zgr_e3} |
---|
362 | Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, |
---|
363 | and (b) analytically derived grid-point position and scale factors. |
---|
364 | For both grids here, the same $w$-point depth has been chosen but in (a) the |
---|
365 | $t$-points are set half way between $w$-points while in (b) they are defined from |
---|
366 | an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$. |
---|
367 | Note the resulting difference between the value of the grid-size $\Delta_k$ and |
---|
368 | those of the scale factor $e_k$. } |
---|
369 | \end{center} \end{figure} |
---|
370 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
371 | |
---|
372 | % ------------------------------------------------------------------------------------------------------------- |
---|
373 | % Choice of horizontal grid |
---|
374 | % ------------------------------------------------------------------------------------------------------------- |
---|
375 | \subsection{Choice of horizontal grid} |
---|
376 | \label{DOM_hgr_msh_choice} |
---|
377 | |
---|
378 | The user has three options available in defining a horizontal grid, which involve |
---|
379 | the namelist variable \np{jphgr\_mesh} of the \ngn{namcfg} namelist. |
---|
380 | \begin{description} |
---|
381 | \item[\np{jphgr\_mesh}=0] The most general curvilinear orthogonal grids. |
---|
382 | The coordinates and their first derivatives with respect to $i$ and $j$ are provided |
---|
383 | in a input file (\ifile{coordinates}), read in \rou{hgr\_read} subroutine of the domhgr module. |
---|
384 | \item[\np{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below). |
---|
385 | For other analytical grids, the \mdl{domhgr} module must be modified by the user. |
---|
386 | \end{description} |
---|
387 | |
---|
388 | There are two simple cases of geographical grids on the sphere. With |
---|
389 | \np{jphgr\_mesh}=1, the grid (expressed in degrees) is regular in space, |
---|
390 | with grid sizes specified by parameters \np{ppe1\_deg} and \np{ppe2\_deg}, |
---|
391 | respectively. Such a geographical grid can be very anisotropic at high latitudes |
---|
392 | because of the convergence of meridians (the zonal scale factors $e_1$ |
---|
393 | become much smaller than the meridional scale factors $e_2$). The Mercator |
---|
394 | grid (\np{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale |
---|
395 | factors in the same way as the zonal ones. In this case, meridional scale factors |
---|
396 | and latitudes are calculated analytically using the formulae appropriate for |
---|
397 | a Mercator projection, based on \np{ppe1\_deg} which is a reference grid spacing |
---|
398 | at the equator (this applies even when the geographical equator is situated outside |
---|
399 | the model domain). |
---|
400 | %%% |
---|
401 | \gmcomment{ give here the analytical expression of the Mercator mesh} |
---|
402 | %%% |
---|
403 | In these two cases (\np{jphgr\_mesh}=1 or 4), the grid position is defined by the |
---|
404 | longitude and latitude of the south-westernmost point (\np{ppglamt0} |
---|
405 | and \np{ppgphi0}). Note that for the Mercator grid the user need only provide |
---|
406 | an approximate starting latitude: the real latitude will be recalculated analytically, |
---|
407 | in order to ensure that the equator corresponds to line passing through $t$- |
---|
408 | and $u$-points. |
---|
409 | |
---|
410 | Rectangular grids ignoring the spherical geometry are defined with |
---|
411 | \np{jphgr\_mesh} = 2, 3, 5. The domain is either an $f$-plane (\np{jphgr\_mesh} = 2, |
---|
412 | Coriolis factor is constant) or a beta-plane (\np{jphgr\_mesh} = 3, the Coriolis factor |
---|
413 | is linear in the $j$-direction). The grid size is uniform in meter in each direction, |
---|
414 | and given by the parameters \np{ppe1\_m} and \np{ppe2\_m} respectively. |
---|
415 | The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero |
---|
416 | with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers, |
---|
417 | and the second $t$-point corresponds to coordinate $gphit=0$. The input |
---|
418 | variable \np{ppglam0} is ignored. \np{ppgphi0} is used to set the reference |
---|
419 | latitude for computation of the Coriolis parameter. In the case of the beta plane, |
---|
420 | \np{ppgphi0} corresponds to the center of the domain. Finally, the special case |
---|
421 | \np{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the |
---|
422 | GYRE configuration, representing a classical mid-latitude double gyre system. |
---|
423 | The rotation allows us to maximize the jet length relative to the gyre areas |
---|
424 | (and the number of grid points). |
---|
425 | |
---|
426 | The choice of the grid must be consistent with the boundary conditions specified |
---|
427 | by \np{jperio}, a parameter found in \ngn{namcfg} namelist (see {\S\ref{LBC}). |
---|
428 | |
---|
429 | % ------------------------------------------------------------------------------------------------------------- |
---|
430 | % Grid files |
---|
431 | % ------------------------------------------------------------------------------------------------------------- |
---|
432 | \subsection{Output Grid files} |
---|
433 | \label{DOM_hgr_files} |
---|
434 | |
---|
435 | All the arrays relating to a particular ocean model configuration (grid-point |
---|
436 | position, scale factors, masks) can be saved in files if $\np{nn\_msh} \not= 0$ |
---|
437 | (namelist variable in \ngn{namdom}). This can be particularly useful for plots and off-line |
---|
438 | diagnostics. In some cases, the user may choose to make a local modification |
---|
439 | of a scale factor in the code. This is the case in global configurations when |
---|
440 | restricting the width of a specific strait (usually a one-grid-point strait that |
---|
441 | happens to be too wide due to insufficient model resolution). An example |
---|
442 | is Gibraltar Strait in the ORCA2 configuration. When such modifications are done, |
---|
443 | the output grid written when $\np{nn\_msh} \not=0$ is no more equal to the input grid. |
---|
444 | |
---|
445 | $\ $\newline % force a new lign |
---|
446 | |
---|
447 | % ================================================================ |
---|
448 | % Domain: Vertical Grid (domzgr) |
---|
449 | % ================================================================ |
---|
450 | \section [Domain: Vertical Grid (\textit{domzgr})] |
---|
451 | {Domain: Vertical Grid \small{(\mdl{domzgr} module)} } |
---|
452 | \label{DOM_zgr} |
---|
453 | %-----------------------------------------nam_zgr & namdom------------------------------------------- |
---|
454 | \namdisplay{namzgr} |
---|
455 | \namdisplay{namdom} |
---|
456 | %------------------------------------------------------------------------------------------------------------- |
---|
457 | |
---|
458 | Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists. |
---|
459 | In the vertical, the model mesh is determined by four things: |
---|
460 | (1) the bathymetry given in meters ; |
---|
461 | (2) the number of levels of the model (\jp{jpk}) ; |
---|
462 | (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors |
---|
463 | (derivatives of the transformation) ; |
---|
464 | and (4) the masking system, $i.e.$ the number of wet model levels at each |
---|
465 | $(i,j)$ column of points. |
---|
466 | |
---|
467 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
468 | \begin{figure}[!tb] \begin{center} |
---|
469 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf} |
---|
470 | \caption{ \label{Fig_z_zps_s_sps} |
---|
471 | The ocean bottom as seen by the model: |
---|
472 | (a) $z$-coordinate with full step, |
---|
473 | (b) $z$-coordinate with partial step, |
---|
474 | (c) $s$-coordinate: terrain following representation, |
---|
475 | (d) hybrid $s-z$ coordinate, |
---|
476 | (e) hybrid $s-z$ coordinate with partial step, and |
---|
477 | (f) same as (e) but with variable volume associated with the non-linear free surface. |
---|
478 | Note that the variable volume option (\key{vvl}) can be used with any of the |
---|
479 | 5 coordinates (a) to (e).} |
---|
480 | \end{center} \end{figure} |
---|
481 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
482 | |
---|
483 | The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters, |
---|
484 | must be done once of all at the beginning of an experiment. It is not intended as an |
---|
485 | option which can be enabled or disabled in the middle of an experiment. Three main |
---|
486 | choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step |
---|
487 | bathymetry (\np{ln\_zco}~=~true), $z$-coordinate with partial step bathymetry |
---|
488 | (\np{ln\_zps}~=~true), or generalized, $s$-coordinate (\np{ln\_sco}~=~true). |
---|
489 | Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate |
---|
490 | (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable |
---|
491 | volume option \key{vvl} ($i.e.$ non-linear free surface), the coordinate follow the |
---|
492 | time-variation of the free surface so that the transformation is time dependent: |
---|
493 | $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step |
---|
494 | bathymetry or $s$-coordinate (hybrid and partial step coordinates have not |
---|
495 | yet been tested in NEMO v2.3). If using $z$-coordinate with partial step bathymetry |
---|
496 | (\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true) |
---|
497 | and partial step are also applied at the ocean/ice shelf interface. |
---|
498 | |
---|
499 | Contrary to the horizontal grid, the vertical grid is computed in the code and no |
---|
500 | provision is made for reading it from a file. The only input file is the bathymetry |
---|
501 | (in meters) (\ifile{bathy\_meter}). |
---|
502 | \footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the |
---|
503 | \ifile{bathy\_meter} file, so that the computation of the number of wet ocean point |
---|
504 | in each water column is by-passed}. |
---|
505 | After reading the bathymetry, the algorithm for vertical grid definition differs |
---|
506 | between the different options: |
---|
507 | \begin{description} |
---|
508 | \item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$. |
---|
509 | \item[\textit{zps}] set a reference coordinate transformation $z_0 (k)$, and |
---|
510 | calculate the thickness of the deepest level at each $(i,j)$ point using the |
---|
511 | bathymetry, to obtain the final three-dimensional depth and scale factor arrays. |
---|
512 | \item[\textit{sco}] smooth the bathymetry to fulfil the hydrostatic consistency |
---|
513 | criteria and set the three-dimensional transformation. |
---|
514 | \item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfil the hydrostatic |
---|
515 | consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and |
---|
516 | possibly introduce masking of extra land points to better fit the original bathymetry file |
---|
517 | \end{description} |
---|
518 | %%% |
---|
519 | \gmcomment{ add the description of the smoothing: envelop topography...} |
---|
520 | %%% |
---|
521 | |
---|
522 | The arrays describing the grid point depths and vertical scale factors |
---|
523 | are three dimensional arrays $(i,j,k)$ even in the case of $z$-coordinate with |
---|
524 | full step bottom topography. In non-linear free surface (\key{vvl}), their knowledge |
---|
525 | is required at \textit{before}, \textit{now} and \textit{after} time step, while they |
---|
526 | do not vary in time in linear free surface case. |
---|
527 | To improve the code readability while providing this flexibility, the vertical coordinate |
---|
528 | and scale factors are defined as functions of |
---|
529 | $(i,j,k)$ with "fs" as prefix (examples: \textit{fse3t\_b, fse3t\_n, fse3t\_a,} |
---|
530 | for the \textit{before}, \textit{now} and \textit{after} scale factors at $t$-point) |
---|
531 | that can be either three different arrays when \key{vvl} is defined, or a single fixed arrays. |
---|
532 | These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory. |
---|
533 | They are used throughout the code, and replaced by the corresponding arrays at |
---|
534 | the time of pre-processing (CPP capability). |
---|
535 | |
---|
536 | % ------------------------------------------------------------------------------------------------------------- |
---|
537 | % Meter Bathymetry |
---|
538 | % ------------------------------------------------------------------------------------------------------------- |
---|
539 | \subsection{Meter Bathymetry} |
---|
540 | \label{DOM_bathy} |
---|
541 | |
---|
542 | Three options are possible for defining the bathymetry, according to the |
---|
543 | namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist): |
---|
544 | \begin{description} |
---|
545 | \item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ |
---|
546 | is given by the coordinate transformation. The domain can either be a closed |
---|
547 | basin or a periodic channel depending on the parameter \np{jperio}. |
---|
548 | \item[\np{nn\_bathy} = -1] a domain with a bump of topography one third of the |
---|
549 | domain width at the central latitude. This is meant for the "EEL-R5" configuration, |
---|
550 | a periodic or open boundary channel with a seamount. |
---|
551 | \item[\np{nn\_bathy} = 1] read a bathymetry and ice shelf draft (if needed). |
---|
552 | The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) |
---|
553 | at each grid point of the model grid. The bathymetry is usually built by interpolating a standard bathymetry product |
---|
554 | ($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also |
---|
555 | defines the coastline: where the bathymetry is zero, no model levels are defined |
---|
556 | (all levels are masked). |
---|
557 | |
---|
558 | The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) |
---|
559 | at each grid point of the model grid. This file is only needed if \np{ln\_isfcav}~=~true. |
---|
560 | Defining the ice shelf draft will also define the ice shelf edge and the grounding line position. |
---|
561 | \end{description} |
---|
562 | |
---|
563 | When a global ocean is coupled to an atmospheric model it is better to represent |
---|
564 | all large water bodies (e.g, great lakes, Caspian sea...) even if the model |
---|
565 | resolution does not allow their communication with the rest of the ocean. |
---|
566 | This is unnecessary when the ocean is forced by fixed atmospheric conditions, |
---|
567 | so these seas can be removed from the ocean domain. The user has the option |
---|
568 | to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the |
---|
569 | code has to be adapted to the user's configuration. |
---|
570 | |
---|
571 | % ------------------------------------------------------------------------------------------------------------- |
---|
572 | % z-coordinate and reference coordinate transformation |
---|
573 | % ------------------------------------------------------------------------------------------------------------- |
---|
574 | \subsection[$z$-coordinate (\np{ln\_zco}] |
---|
575 | {$z$-coordinate (\np{ln\_zco}=true) and reference coordinate} |
---|
576 | \label{DOM_zco} |
---|
577 | |
---|
578 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
579 | \begin{figure}[!tb] \begin{center} |
---|
580 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf} |
---|
581 | \caption{ \label{Fig_zgr} |
---|
582 | Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for |
---|
583 | (a) T-point depth and (b) the associated scale factor as computed |
---|
584 | from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.} |
---|
585 | \end{center} \end{figure} |
---|
586 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
587 | |
---|
588 | The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ |
---|
589 | and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on |
---|
590 | Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the |
---|
591 | ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the |
---|
592 | additional $t$-point at $jk=jpk$ is below the sea floor and is not used. |
---|
593 | The vertical location of $w$- and $t$-levels is defined from the analytic expression |
---|
594 | of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the |
---|
595 | vertical scale factors. The user must provide the analytical expression of both |
---|
596 | $z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} |
---|
597 | through statement functions, using parameters provided in the \ngn{namcfg} namelist. |
---|
598 | |
---|
599 | It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr=0}). |
---|
600 | In that case, the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} |
---|
601 | (total ocean depth in meters) fully define the grid. |
---|
602 | |
---|
603 | For climate-related studies it is often desirable to concentrate the vertical resolution |
---|
604 | near the ocean surface. The following function is proposed as a standard for a |
---|
605 | $z$-coordinate (with either full or partial steps): |
---|
606 | \begin{equation} \label{DOM_zgr_ana} |
---|
607 | \begin{split} |
---|
608 | z_0 (k) &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ |
---|
609 | e_3^0 (k) &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right| |
---|
610 | \end{split} |
---|
611 | \end{equation} |
---|
612 | where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an |
---|
613 | expression allows us to define a nearly uniform vertical location of levels at the |
---|
614 | ocean top and bottom with a smooth hyperbolic tangent transition in between |
---|
615 | (Fig.~\ref{Fig_zgr}). |
---|
616 | |
---|
617 | If the ice shelf cavities are opened (\np{ln\_isfcav}=~true~}), the definition of $z_0$ is the same. |
---|
618 | However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to: |
---|
619 | \begin{equation} \label{DOM_zgr_ana} |
---|
620 | \begin{split} |
---|
621 | e_3^T(k) &= z_W (k+1) - z_W (k) \\ |
---|
622 | e_3^W(k) &= z_T (k) - z_T (k-1) \\ |
---|
623 | \end{split} |
---|
624 | \end{equation} |
---|
625 | This formulation decrease the self-generated circulation into the ice shelf cavity |
---|
626 | (which can, in extreme case, leads to blow up).\\ |
---|
627 | |
---|
628 | |
---|
629 | The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the |
---|
630 | surface (bottom) layers and a depth which varies from 0 at the sea surface to a |
---|
631 | minimum of $-5000~m$. This leads to the following conditions: |
---|
632 | \begin{equation} \label{DOM_zgr_coef} |
---|
633 | \begin{split} |
---|
634 | e_3 (1+1/2) &=10. \\ |
---|
635 | e_3 (jpk-1/2) &=500. \\ |
---|
636 | z(1) &=0. \\ |
---|
637 | z(jpk) &=-5000. \\ |
---|
638 | \end{split} |
---|
639 | \end{equation} |
---|
640 | |
---|
641 | With the choice of the stretching $h_{cr} =3$ and the number of levels |
---|
642 | \jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in |
---|
643 | \eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is |
---|
644 | satisfied, through an optimisation procedure using a bisection method. For the first |
---|
645 | standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$, |
---|
646 | $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and |
---|
647 | scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and |
---|
648 | given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters |
---|
649 | \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist. |
---|
650 | |
---|
651 | Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is |
---|
652 | possible to recalculate them. In that case the user sets |
---|
653 | \np{ppsur}=\np{ppa0}=\np{ppa1}=999999., in \ngn{namcfg} namelist, |
---|
654 | and specifies instead the four following parameters: |
---|
655 | \begin{itemize} |
---|
656 | \item \np{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger |
---|
657 | \np{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual. |
---|
658 | \item \np{ppkth}=$h_{th} $: is approximately the model level at which maximum |
---|
659 | stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) |
---|
660 | \item \np{ppdzmin}: minimum thickness for the top layer (in meters) |
---|
661 | \item \np{pphmax}: total depth of the ocean (meters). |
---|
662 | \end{itemize} |
---|
663 | As an example, for the $45$ layers used in the DRAKKAR configuration those |
---|
664 | parameters are: \jp{jpk}=46, \np{ppacr}=9, \np{ppkth}=23.563, \np{ppdzmin}=6m, |
---|
665 | \np{pphmax}=5750m. |
---|
666 | |
---|
667 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
668 | \begin{table} \begin{center} \begin{tabular}{c||r|r|r|r} |
---|
669 | \hline |
---|
670 | \textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w } \\ \hline |
---|
671 | 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
672 | 2 & \textbf{15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
673 | 3 & \textbf{25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
674 | 4 & \textbf{35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ \hline |
---|
675 | 5 & \textbf{45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ \hline |
---|
676 | 6 & \textbf{55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ \hline |
---|
677 | 7 & \textbf{65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ \hline |
---|
678 | 8 & \textbf{75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ \hline |
---|
679 | 9 & \textbf{85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ \hline |
---|
680 | 10 & \textbf{95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ \hline |
---|
681 | 11 & \textbf{105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ \hline |
---|
682 | 12 & \textbf{116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ \hline |
---|
683 | 13 & \textbf{128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ \hline |
---|
684 | 14 & \textbf{142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ \hline |
---|
685 | 15 & \textbf{158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ \hline |
---|
686 | 16 & \textbf{181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ \hline |
---|
687 | 17 & \textbf{216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ \hline |
---|
688 | 18 & \textbf{272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ \hline |
---|
689 | 19 & \textbf{364.30} & 312.74 & \textbf{116.11} & 90.99 \\ \hline |
---|
690 | 20 & \textbf{511.53} & 429.72 & \textbf{181.55} & 146.43 \\ \hline |
---|
691 | 21 & \textbf{732.20} & 611.89 & \textbf{261.03} & 220.35 \\ \hline |
---|
692 | 22 & \textbf{1033.22}& 872.87 & \textbf{339.39} & 301.42 \\ \hline |
---|
693 | 23 & \textbf{1405.70}& 1211.59 & \textbf{402.26} & 373.31 \\ \hline |
---|
694 | 24 & \textbf{1830.89}& 1612.98 & \textbf{444.87} & 426.00 \\ \hline |
---|
695 | 25 & \textbf{2289.77}& 2057.13 & \textbf{470.55} & 459.47 \\ \hline |
---|
696 | 26 & \textbf{2768.24}& 2527.22 & \textbf{484.95} & 478.83 \\ \hline |
---|
697 | 27 & \textbf{3257.48}& 3011.90 & \textbf{492.70} & 489.44 \\ \hline |
---|
698 | 28 & \textbf{3752.44}& 3504.46 & \textbf{496.78} & 495.07 \\ \hline |
---|
699 | 29 & \textbf{4250.40}& 4001.16 & \textbf{498.90} & 498.02 \\ \hline |
---|
700 | 30 & \textbf{4749.91}& 4500.02 & \textbf{500.00} & 499.54 \\ \hline |
---|
701 | 31 & \textbf{5250.23}& 5000.00 & \textbf{500.56} & 500.33 \\ \hline |
---|
702 | \end{tabular} \end{center} |
---|
703 | \caption{ \label{Tab_orca_zgr} |
---|
704 | Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed |
---|
705 | from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}} |
---|
706 | \end{table} |
---|
707 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
708 | |
---|
709 | % ------------------------------------------------------------------------------------------------------------- |
---|
710 | % z-coordinate with partial step |
---|
711 | % ------------------------------------------------------------------------------------------------------------- |
---|
712 | \subsection [$z$-coordinate with partial step (\np{ln\_zps})] |
---|
713 | {$z$-coordinate with partial step (\np{ln\_zps}=.true.)} |
---|
714 | \label{DOM_zps} |
---|
715 | %--------------------------------------------namdom------------------------------------------------------- |
---|
716 | \namdisplay{namdom} |
---|
717 | %-------------------------------------------------------------------------------------------------------------- |
---|
718 | |
---|
719 | In $z$-coordinate partial step, the depths of the model levels are defined by the |
---|
720 | reference analytical function $z_0 (k)$ as described in the previous |
---|
721 | section, \emph{except} in the bottom layer. The thickness of the bottom layer is |
---|
722 | allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a |
---|
723 | better representation of the bathymetry, especially in the case of small |
---|
724 | slopes (where the bathymetry varies by less than one level thickness from |
---|
725 | one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been |
---|
726 | defined in the absence of bathymetry. With partial steps, layers from 1 to |
---|
727 | \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1) |
---|
728 | is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the |
---|
729 | maximum thickness allowed is $2*e_{3t}(jpk-1)$. This has to be kept in mind when |
---|
730 | specifying values in \ngn{namdom} namelist, as the maximum depth \np{pphmax} |
---|
731 | in partial steps: for example, with |
---|
732 | \np{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth |
---|
733 | allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). |
---|
734 | Two variables in the namdom namelist are used to define the partial step |
---|
735 | vertical grid. The mimimum water thickness (in meters) allowed for a cell |
---|
736 | partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min} |
---|
737 | (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction, |
---|
738 | usually 10\%, of the default thickness $e_{3t}(jk)$). |
---|
739 | |
---|
740 | \gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } } |
---|
741 | |
---|
742 | % ------------------------------------------------------------------------------------------------------------- |
---|
743 | % s-coordinate |
---|
744 | % ------------------------------------------------------------------------------------------------------------- |
---|
745 | \subsection [$s$-coordinate (\np{ln\_sco})] |
---|
746 | {$s$-coordinate (\np{ln\_sco}=true)} |
---|
747 | \label{DOM_sco} |
---|
748 | %------------------------------------------nam_zgr_sco--------------------------------------------------- |
---|
749 | \namdisplay{namzgr_sco} |
---|
750 | %-------------------------------------------------------------------------------------------------------------- |
---|
751 | Options are defined in \ngn{namzgr\_sco}. |
---|
752 | In $s$-coordinate (\np{ln\_sco}~=~true), the depth and thickness of the model |
---|
753 | levels are defined from the product of a depth field and either a stretching |
---|
754 | function or its derivative, respectively: |
---|
755 | |
---|
756 | \begin{equation} \label{DOM_sco_ana} |
---|
757 | \begin{split} |
---|
758 | z(k) &= h(i,j) \; z_0(k) \\ |
---|
759 | e_3(k) &= h(i,j) \; z_0'(k) |
---|
760 | \end{split} |
---|
761 | \end{equation} |
---|
762 | |
---|
763 | where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point |
---|
764 | location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea |
---|
765 | surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean |
---|
766 | depth, since a mixed step-like and bottom-following representation of the |
---|
767 | topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f). |
---|
768 | The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. |
---|
769 | |
---|
770 | Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application. |
---|
771 | |
---|
772 | The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true (\np{ln\_sco\_SH94}~=~false and \np{ln\_sco\_SF12}~=~false.) This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}: |
---|
773 | |
---|
774 | \begin{equation} |
---|
775 | z = s_{min}+C\left(s\right)\left(H-s_{min}\right) |
---|
776 | \label{eq:SH94_1} |
---|
777 | \end{equation} |
---|
778 | |
---|
779 | where $s_{min}$ is the depth at which the s-coordinate stretching starts and allows a z-coordinate to placed on top of the stretched coordinate, and z is the depth (negative down from the asea surface). |
---|
780 | |
---|
781 | \begin{equation} |
---|
782 | s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1 |
---|
783 | \label{eq:s} |
---|
784 | \end{equation} |
---|
785 | |
---|
786 | \begin{equation} \label{DOM_sco_function} |
---|
787 | \begin{split} |
---|
788 | C(s) &= \frac{ \left[ \tanh{ \left( \theta \, (s+b) \right)} |
---|
789 | - \tanh{ \left( \theta \, b \right)} \right]} |
---|
790 | {2\;\sinh \left( \theta \right)} |
---|
791 | \end{split} |
---|
792 | \end{equation} |
---|
793 | |
---|
794 | A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling: |
---|
795 | |
---|
796 | \begin{equation} |
---|
797 | C\left(s\right) = \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} + \\ |
---|
798 | b\frac{ \tanh \left[ \theta \left(s + \frac{1}{2} \right)\right] - \tanh\left(\frac{\theta}{2}\right)}{ 2\tanh\left (\frac{\theta}{2}\right)} |
---|
799 | \label{eq:SH94_2} |
---|
800 | \end{equation} |
---|
801 | |
---|
802 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
803 | \begin{figure}[!ht] \begin{center} |
---|
804 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf} |
---|
805 | \caption{ \label{Fig_sco_function} |
---|
806 | Examples of the stretching function applied to a seamount; from left to right: |
---|
807 | surface, surface and bottom, and bottom intensified resolutions} |
---|
808 | \end{center} \end{figure} |
---|
809 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
810 | |
---|
811 | where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and |
---|
812 | bottom control parameters such that $0\leqslant \theta \leqslant 20$, and |
---|
813 | $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom |
---|
814 | increase of the vertical resolution (Fig.~\ref{Fig_sco_function}). |
---|
815 | |
---|
816 | Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that: |
---|
817 | |
---|
818 | \begin{equation} |
---|
819 | z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1 |
---|
820 | \label{eq:z} |
---|
821 | \end{equation} |
---|
822 | |
---|
823 | The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate: |
---|
824 | |
---|
825 | \begin{equation} \label{DOM_gamma_deriv} |
---|
826 | \gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right) |
---|
827 | \end{equation} |
---|
828 | |
---|
829 | Where: |
---|
830 | \begin{equation} \label{DOM_gamma} |
---|
831 | f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2} \quad \text{ and } \quad \sigma = \frac{k}{n-1} |
---|
832 | \end{equation} |
---|
833 | |
---|
834 | This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth: |
---|
835 | |
---|
836 | \begin{equation} \label{DOM_zb} |
---|
837 | Z_b= h a + b |
---|
838 | \end{equation} |
---|
839 | |
---|
840 | where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively. |
---|
841 | |
---|
842 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
843 | \begin{figure}[!ht] |
---|
844 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/FIG_DOM_compare_coordinates_surface.pdf} |
---|
845 | \caption{A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.} |
---|
846 | \label{fig_compare_coordinates_surface} |
---|
847 | \end{figure} |
---|
848 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
849 | |
---|
850 | This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell thicknesses in real space. This is not to be confused with the hybrid schemes that superimpose geopotential coordinates on terrain following coordinates thus creating a non-analytical vertical coordinate that therefore may suffer from large gradients in the vertical resolutions. This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes. |
---|
851 | |
---|
852 | As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than the critical depth $h_c$. In this example two options are available in depths shallower than $h_c$, with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$. |
---|
853 | |
---|
854 | Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as large slopes lead to hydrostatic consistency. A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, and is output as part of the model mesh file at the start of the run. |
---|
855 | |
---|
856 | % ------------------------------------------------------------------------------------------------------------- |
---|
857 | % z*- or s*-coordinate |
---|
858 | % ------------------------------------------------------------------------------------------------------------- |
---|
859 | \subsection{$z^*$- or $s^*$-coordinate (add \key{vvl}) } |
---|
860 | \label{DOM_zgr_vvl} |
---|
861 | |
---|
862 | This option is described in the Report by Levier \textit{et al.} (2007), available on |
---|
863 | the \NEMO web site. |
---|
864 | |
---|
865 | %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances |
---|
866 | |
---|
867 | % ------------------------------------------------------------------------------------------------------------- |
---|
868 | % level bathymetry and mask |
---|
869 | % ------------------------------------------------------------------------------------------------------------- |
---|
870 | \subsection{level bathymetry and mask} |
---|
871 | \label{DOM_msk} |
---|
872 | |
---|
873 | Whatever the vertical coordinate used, the model offers the possibility of |
---|
874 | representing the bottom topography with steps that follow the face of the |
---|
875 | model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of |
---|
876 | the steps in the horizontal is defined in a 2D integer array, mbathy, which |
---|
877 | gives the number of ocean levels ($i.e.$ those that are not masked) at each |
---|
878 | $t$-point. mbathy is computed from the meter bathymetry using the definiton of |
---|
879 | gdept as the number of $t$-points which gdept $\leq$ bathy. |
---|
880 | |
---|
881 | Modifications of the model bathymetry are performed in the \textit{bat\_ctl} |
---|
882 | routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points |
---|
883 | that do not communicate with another ocean point at the same level are eliminated.\\ |
---|
884 | |
---|
885 | As for the representation of bathymetry, a 2D integer array, misfdep, is created. |
---|
886 | misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. |
---|
887 | By default, misfdep(:,:)=1 and no cells are masked. |
---|
888 | |
---|
889 | In case of ice shelf cavities (\np{ln\_isfcav}~=~true), modifications of the model bathymetry and ice shelf draft in |
---|
890 | the cavities are performed through the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked: |
---|
891 | if only one cell on the water column is opened at $t$-, $u$- or $v$-points, the bathymetry or the ice shelf draft is dug to have a 2-level water column |
---|
892 | (i.e. two unmasked levels). If the incompatibility is too strong (i.e. need to dig more than one cell), the entire water column is masked.\\ |
---|
893 | |
---|
894 | From the \textit{mbathy} array, the mask fields are defined as follows: |
---|
895 | \begin{align*} |
---|
896 | tmask(i,j,k) &= \begin{cases} \; 0& \text{ if $k < misfdep(i,j) $ } \\ |
---|
897 | \; 1& \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$ } \\ |
---|
898 | \; 0& \text{ if $k > mbathy(i,j)$ } \end{cases} \\ |
---|
899 | umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ |
---|
900 | vmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\ |
---|
901 | fmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ |
---|
902 | & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ |
---|
903 | wmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1) |
---|
904 | \end{align*} |
---|
905 | |
---|
906 | Note, wmask is now defined. It allows, in case of ice shelves, |
---|
907 | to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. |
---|
908 | |
---|
909 | The specification of closed lateral boundaries requires that at least the first and last |
---|
910 | rows and columns of the \textit{mbathy} array are set to zero. In the particular |
---|
911 | case of an east-west cyclical boundary condition, \textit{mbathy} has its last |
---|
912 | column equal to the second one and its first column equal to the last but one |
---|
913 | (and so too the mask arrays) (see \S~\ref{LBC_jperio}). |
---|
914 | |
---|
915 | %%% |
---|
916 | \gmcomment{ \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}. } |
---|
917 | %%% |
---|
918 | |
---|
919 | % ================================================================ |
---|
920 | % Domain: Initial State (dtatsd & istate) |
---|
921 | % ================================================================ |
---|
922 | \section [Domain: Initial State (\textit{istate and dtatsd})] |
---|
923 | {Domain: Initial State \small{(\mdl{istate} and \mdl{dtatsd} modules)} } |
---|
924 | \label{DTA_tsd} |
---|
925 | %-----------------------------------------namtsd------------------------------------------- |
---|
926 | \namdisplay{namtsd} |
---|
927 | %------------------------------------------------------------------------------------------ |
---|
928 | |
---|
929 | Options are defined in \ngn{namtsd}. |
---|
930 | By default, the ocean start from rest (the velocity field is set to zero) and the initialization of |
---|
931 | temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter. |
---|
932 | \begin{description} |
---|
933 | \item[ln\_tsd\_init = .true.] use a T and S input files that can be given on the model grid itself or |
---|
934 | on their native input data grid. In the latter case, the data will be interpolated on-the-fly both in the |
---|
935 | horizontal and the vertical to the model grid (see \S~\ref{SBC_iof}). The information relative to the |
---|
936 | input files are given in the \np{sn\_tem} and \np{sn\_sal} structures. |
---|
937 | The computation is done in the \mdl{dtatsd} module. |
---|
938 | \item[ln\_tsd\_init = .false.] use constant salinity value of 35.5 psu and an analytical profile of temperature |
---|
939 | (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module. |
---|
940 | \end{description} |
---|