1 | % ================================================================ |
---|
2 | % Chapter Ñ Lateral Boundary Condition (LBC) |
---|
3 | % ================================================================ |
---|
4 | \chapter{Lateral Boundary Condition (LBC) } |
---|
5 | \label{LBC} |
---|
6 | \minitoc |
---|
7 | |
---|
8 | \newpage |
---|
9 | $\ $\newline % force a new ligne |
---|
10 | |
---|
11 | |
---|
12 | %gm% add here introduction to this chapter |
---|
13 | |
---|
14 | % ================================================================ |
---|
15 | % Boundary Condition at the Coast |
---|
16 | % ================================================================ |
---|
17 | \section{Boundary Condition at the Coast (\np{rn\_shlat})} |
---|
18 | \label{LBC_coast} |
---|
19 | %--------------------------------------------nam_lbc------------------------------------------------------- |
---|
20 | \namdisplay{namlbc} |
---|
21 | %-------------------------------------------------------------------------------------------------------------- |
---|
22 | |
---|
23 | %The lateral ocean boundary conditions contiguous to coastlines are Neumann conditions for heat and salt (no flux across boundaries) and Dirichlet conditions for momentum (ranging from free-slip to "strong" no-slip). They are handled automatically by the mask system (see \S\ref{DOM_msk}). |
---|
24 | |
---|
25 | %OPA allows land and topography grid points in the computational domain due to the presence of continents or islands, and includes the use of a full or partial step representation of bottom topography. The computation is performed over the whole domain, i.e. we do not try to restrict the computation to ocean-only points. This choice has two motivations. Firstly, working on ocean only grid points overloads the code and harms the code readability. Secondly, and more importantly, it drastically reduces the vector portion of the computation, leading to a dramatic increase of CPU time requirement on vector computers. The current section describes how the masking affects the computation of the various terms of the equations with respect to the boundary condition at solid walls. The process of defining which areas are to be masked is described in \S\ref{DOM_msk}. |
---|
26 | |
---|
27 | The discrete representation of a domain with complex boundaries (coastlines and |
---|
28 | bottom topography) leads to arrays that include large portions where a computation |
---|
29 | is not required as the model variables remain at zero. Nevertheless, vectorial |
---|
30 | supercomputers are far more efficient when computing over a whole array, and the |
---|
31 | readability of a code is greatly improved when boundary conditions are applied in |
---|
32 | an automatic way rather than by a specific computation before or after each |
---|
33 | computational loop. An efficient way to work over the whole domain while specifying |
---|
34 | the boundary conditions, is to use multiplication by mask arrays in the computation. |
---|
35 | A mask array is a matrix whose elements are $1$ in the ocean domain and $0$ |
---|
36 | elsewhere. A simple multiplication of a variable by its own mask ensures that it will |
---|
37 | remain zero over land areas. Since most of the boundary conditions consist of a |
---|
38 | zero flux across the solid boundaries, they can be simply applied by multiplying |
---|
39 | variables by the correct mask arrays, $i.e.$ the mask array of the grid point where |
---|
40 | the flux is evaluated. For example, the heat flux in the \textbf{i}-direction is evaluated |
---|
41 | at $u$-points. Evaluating this quantity as, |
---|
42 | |
---|
43 | \begin{equation} \label{Eq_lbc_aaaa} |
---|
44 | \frac{A^{lT} }{e_1 }\frac{\partial T}{\partial i}\equiv \frac{A_u^{lT} |
---|
45 | }{e_{1u} } \; \delta _{i+1 / 2} \left[ T \right]\;\;mask_u |
---|
46 | \end{equation} |
---|
47 | (where mask$_{u}$ is the mask array at a $u$-point) ensures that the heat flux is |
---|
48 | zero inside land and at the boundaries, since mask$_{u}$ is zero at solid boundaries |
---|
49 | which in this case are defined at $u$-points (normal velocity $u$ remains zero at |
---|
50 | the coast) (Fig.~\ref{Fig_LBC_uv}). |
---|
51 | |
---|
52 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
53 | \begin{figure}[!t] \begin{center} |
---|
54 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_uv.pdf} |
---|
55 | \caption{ \label{Fig_LBC_uv} |
---|
56 | Lateral boundary (thick line) at T-level. The velocity normal to the boundary is set to zero.} |
---|
57 | \end{center} \end{figure} |
---|
58 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
59 | |
---|
60 | For momentum the situation is a bit more complex as two boundary conditions |
---|
61 | must be provided along the coast (one each for the normal and tangential velocities). |
---|
62 | The boundary of the ocean in the C-grid is defined by the velocity-faces. |
---|
63 | For example, at a given $T$-level, the lateral boundary (a coastline or an intersection |
---|
64 | with the bottom topography) is made of segments joining $f$-points, and normal |
---|
65 | velocity points are located between two $f-$points (Fig.~\ref{Fig_LBC_uv}). |
---|
66 | The boundary condition on the normal velocity (no flux through solid boundaries) |
---|
67 | can thus be easily implemented using the mask system. The boundary condition |
---|
68 | on the tangential velocity requires a more specific treatment. This boundary |
---|
69 | condition influences the relative vorticity and momentum diffusive trends, and is |
---|
70 | required in order to compute the vorticity at the coast. Four different types of |
---|
71 | lateral boundary condition are available, controlled by the value of the \np{rn\_shlat} |
---|
72 | namelist parameter. (The value of the mask$_{f}$ array along the coastline is set |
---|
73 | equal to this parameter.) These are: |
---|
74 | |
---|
75 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
76 | \begin{figure}[!p] \begin{center} |
---|
77 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_shlat.pdf} |
---|
78 | \caption{ \label{Fig_LBC_shlat} |
---|
79 | lateral boundary condition (a) free-slip ($rn\_shlat=0$) ; (b) no-slip ($rn\_shlat=2$) |
---|
80 | ; (c) "partial" free-slip ($0<rn\_shlat<2$) and (d) "strong" no-slip ($2<rn\_shlat$). |
---|
81 | Implied "ghost" velocity inside land area is display in grey. } |
---|
82 | \end{center} \end{figure} |
---|
83 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
84 | |
---|
85 | \begin{description} |
---|
86 | |
---|
87 | \item[free-slip boundary condition (\np{rn\_shlat}=0): ] the tangential velocity at the |
---|
88 | coastline is equal to the offshore velocity, $i.e.$ the normal derivative of the |
---|
89 | tangential velocity is zero at the coast, so the vorticity: mask$_{f}$ array is set |
---|
90 | to zero inside the land and just at the coast (Fig.~\ref{Fig_LBC_shlat}-a). |
---|
91 | |
---|
92 | \item[no-slip boundary condition (\np{rn\_shlat}=2): ] the tangential velocity vanishes |
---|
93 | at the coastline. Assuming that the tangential velocity decreases linearly from |
---|
94 | the closest ocean velocity grid point to the coastline, the normal derivative is |
---|
95 | evaluated as if the velocities at the closest land velocity gridpoint and the closest |
---|
96 | ocean velocity gridpoint were of the same magnitude but in the opposite direction |
---|
97 | (Fig.~\ref{Fig_LBC_shlat}-b). Therefore, the vorticity along the coastlines is given by: |
---|
98 | |
---|
99 | \begin{equation*} |
---|
100 | \zeta \equiv 2 \left(\delta_{i+1/2} \left[e_{2v} v \right] - \delta_{j+1/2} \left[e_{1u} u \right] \right) / \left(e_{1f} e_{2f} \right) \ , |
---|
101 | \end{equation*} |
---|
102 | where $u$ and $v$ are masked fields. Setting the mask$_{f}$ array to $2$ along |
---|
103 | the coastline provides a vorticity field computed with the no-slip boundary condition, |
---|
104 | simply by multiplying it by the mask$_{f}$ : |
---|
105 | \begin{equation} \label{Eq_lbc_bbbb} |
---|
106 | \zeta \equiv \frac{1}{e_{1f} {\kern 1pt}e_{2f} }\left( {\delta _{i+1/2} |
---|
107 | \left[ {e_{2v} \,v} \right]-\delta _{j+1/2} \left[ {e_{1u} \,u} \right]} |
---|
108 | \right)\;\mbox{mask}_f |
---|
109 | \end{equation} |
---|
110 | |
---|
111 | \item["partial" free-slip boundary condition (0$<$\np{rn\_shlat}$<$2): ] the tangential |
---|
112 | velocity at the coastline is smaller than the offshore velocity, $i.e.$ there is a lateral |
---|
113 | friction but not strong enough to make the tangential velocity at the coast vanish |
---|
114 | (Fig.~\ref{Fig_LBC_shlat}-c). This can be selected by providing a value of mask$_{f}$ |
---|
115 | strictly inbetween $0$ and $2$. |
---|
116 | |
---|
117 | \item["strong" no-slip boundary condition (2$<$\np{rn\_shlat}): ] the viscous boundary |
---|
118 | layer is assumed to be smaller than half the grid size (Fig.~\ref{Fig_LBC_shlat}-d). |
---|
119 | The friction is thus larger than in the no-slip case. |
---|
120 | |
---|
121 | \end{description} |
---|
122 | |
---|
123 | Note that when the bottom topography is entirely represented by the $s$-coor-dinates |
---|
124 | (pure $s$-coordinate), the lateral boundary condition on tangential velocity is of much |
---|
125 | less importance as it is only applied next to the coast where the minimum water depth |
---|
126 | can be quite shallow. |
---|
127 | |
---|
128 | The alternative numerical implementation of the no-slip boundary conditions for an |
---|
129 | arbitrary coast line of \citet{Shchepetkin1996} is also available through the |
---|
130 | \key{noslip\_accurate} CPP key. It is based on a fourth order evaluation of the shear at the |
---|
131 | coast which, in turn, allows a true second order scheme in the interior of the domain |
---|
132 | ($i.e.$ the numerical boundary scheme simulates the truncation error of the numerical |
---|
133 | scheme used in the interior of the domain). \citet{Shchepetkin1996} found that such a |
---|
134 | technique considerably improves the quality of the numerical solution. In \NEMO, such |
---|
135 | spectacular improvements have not been found in the half-degree global ocean |
---|
136 | (ORCA05), but significant reductions of numerically induced coastal upwellings were |
---|
137 | found in an eddy resolving simulation of the Alboran Sea \citep{Olivier_PhD01}. |
---|
138 | Nevertheless, since a no-slip boundary condition is not recommended in an eddy |
---|
139 | permitting or resolving simulation \citep{Penduff_al_OS07}, the use of this option is also |
---|
140 | not recommended. |
---|
141 | |
---|
142 | In practice, the no-slip accurate option changes the way the curl is evaluated at the |
---|
143 | coast (see \mdl{divcur} module), and requires the nature of each coastline grid point |
---|
144 | (convex or concave corners, straight north-south or east-west coast) to be specified. |
---|
145 | This is performed in routine \rou{dom\_msk\_nsa} in the \mdl{domask} module. |
---|
146 | |
---|
147 | % ================================================================ |
---|
148 | % Boundary Condition around the Model Domain |
---|
149 | % ================================================================ |
---|
150 | \section{Model Domain Boundary Condition (\jp{jperio})} |
---|
151 | \label{LBC_jperio} |
---|
152 | |
---|
153 | At the model domain boundaries several choices are offered: closed, cyclic east-west, |
---|
154 | south symmetric across the equator, a north-fold, and combination closed-north fold |
---|
155 | or cyclic-north-fold. The north-fold boundary condition is associated with the 3-pole ORCA mesh. |
---|
156 | |
---|
157 | % ------------------------------------------------------------------------------------------------------------- |
---|
158 | % Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2) |
---|
159 | % ------------------------------------------------------------------------------------------------------------- |
---|
160 | \subsection{Closed, cyclic, south symmetric (\jp{jperio} = 0, 1 or 2)} |
---|
161 | \label{LBC_jperio012} |
---|
162 | |
---|
163 | The choice of closed, cyclic or symmetric model domain boundary condition is made |
---|
164 | by setting \jp{jperio} to 0, 1 or 2 in file \mdl{par\_oce}. Each time such a boundary |
---|
165 | condition is needed, it is set by a call to routine \mdl{lbclnk}. The computation of |
---|
166 | momentum and tracer trends proceeds from $i=2$ to $i=jpi-1$ and from $j=2$ to |
---|
167 | $j=jpj-1$, $i.e.$ in the model interior. To choose a lateral model boundary condition |
---|
168 | is to specify the first and last rows and columns of the model variables. |
---|
169 | |
---|
170 | \begin{description} |
---|
171 | |
---|
172 | \item[For closed boundary (\textit{jperio=0})], solid walls are imposed at all model |
---|
173 | boundaries: first and last rows and columns are set to zero. |
---|
174 | |
---|
175 | \item[For cyclic east-west boundary (\textit{jperio=1})], first and last rows are set |
---|
176 | to zero (closed) whilst the first column is set to the value of the last-but-one column |
---|
177 | and the last column to the value of the second one (Fig.~\ref{Fig_LBC_jperio}-a). |
---|
178 | Whatever flows out of the eastern (western) end of the basin enters the western |
---|
179 | (eastern) end. Note that there is no option for north-south cyclic or for doubly |
---|
180 | cyclic cases. |
---|
181 | |
---|
182 | \item[For symmetric boundary condition across the equator (\textit{jperio=2})], |
---|
183 | last rows, and first and last columns are set to zero (closed). The row of symmetry |
---|
184 | is chosen to be the $u$- and $T-$points equator line ($j=2$, i.e. at the southern |
---|
185 | end of the domain). For arrays defined at $u-$ or $T-$points, the first row is set |
---|
186 | to the value of the third row while for most of $v$- and $f$-point arrays ($v$, $\zeta$, |
---|
187 | $j\psi$, but \gmcomment{not sure why this is "but"} scalar arrays such as eddy coefficients) |
---|
188 | the first row is set to minus the value of the second row (Fig.~\ref{Fig_LBC_jperio}-b). |
---|
189 | Note that this boundary condition is not yet available for the case of a massively |
---|
190 | parallel computer (\textbf{key{\_}mpp} defined). |
---|
191 | |
---|
192 | \end{description} |
---|
193 | |
---|
194 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
195 | \begin{figure}[!t] \begin{center} |
---|
196 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_jperio.pdf} |
---|
197 | \caption{ \label{Fig_LBC_jperio} |
---|
198 | setting of (a) east-west cyclic (b) symmetric across the equator boundary conditions.} |
---|
199 | \end{center} \end{figure} |
---|
200 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
201 | |
---|
202 | % ------------------------------------------------------------------------------------------------------------- |
---|
203 | % North fold (\textit{jperio = 3 }to $6)$ |
---|
204 | % ------------------------------------------------------------------------------------------------------------- |
---|
205 | \subsection{North-fold (\textit{jperio = 3 }to $6)$ } |
---|
206 | \label{LBC_north_fold} |
---|
207 | |
---|
208 | The north fold boundary condition has been introduced in order to handle the north |
---|
209 | boundary of a three-polar ORCA grid. Such a grid has two poles in the northern hemisphere. |
---|
210 | \colorbox{yellow}{to be completed...} |
---|
211 | |
---|
212 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
213 | \begin{figure}[!t] \begin{center} |
---|
214 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_North_Fold_T.pdf} |
---|
215 | \caption{ \label{Fig_North_Fold_T} |
---|
216 | North fold boundary with a $T$-point pivot and cyclic east-west boundary condition |
---|
217 | ($jperio=4$), as used in ORCA 2, 1/4, and 1/12. Pink shaded area corresponds |
---|
218 | to the inner domain mask (see text). } |
---|
219 | \end{center} \end{figure} |
---|
220 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
221 | |
---|
222 | % ==================================================================== |
---|
223 | % Exchange with neighbouring processors |
---|
224 | % ==================================================================== |
---|
225 | \section [Exchange with neighbouring processors (\textit{lbclnk}, \textit{lib\_mpp})] |
---|
226 | {Exchange with neighbouring processors (\mdl{lbclnk}, \mdl{lib\_mpp})} |
---|
227 | \label{LBC_mpp} |
---|
228 | |
---|
229 | For massively parallel processing (mpp), a domain decomposition method is used. |
---|
230 | The basic idea of the method is to split the large computation domain of a numerical |
---|
231 | experiment into several smaller domains and solve the set of equations by addressing |
---|
232 | independent local problems. Each processor has its own local memory and computes |
---|
233 | the model equation over a subdomain of the whole model domain. The subdomain |
---|
234 | boundary conditions are specified through communications between processors |
---|
235 | which are organized by explicit statements (message passing method). |
---|
236 | |
---|
237 | A big advantage is that the method does not need many modifications of the initial |
---|
238 | FORTRAN code. From the modeller's point of view, each sub domain running on |
---|
239 | a processor is identical to the "mono-domain" code. In addition, the programmer |
---|
240 | manages the communications between subdomains, and the code is faster when |
---|
241 | the number of processors is increased. The porting of OPA code on an iPSC860 |
---|
242 | was achieved during Guyon's PhD [Guyon et al. 1994, 1995] in collaboration with |
---|
243 | CETIIS and ONERA. The implementation in the operational context and the studies |
---|
244 | of performance on a T3D and T3E Cray computers have been made in collaboration |
---|
245 | with IDRIS and CNRS. The present implementation is largely inspired by Guyon's |
---|
246 | work [Guyon 1995]. |
---|
247 | |
---|
248 | The parallelization strategy is defined by the physical characteristics of the |
---|
249 | ocean model. Second order finite difference schemes lead to local discrete |
---|
250 | operators that depend at the very most on one neighbouring point. The only |
---|
251 | non-local computations concern the vertical physics (implicit diffusion, 1.5 |
---|
252 | turbulent closure scheme, ...) (delocalization over the whole water column), |
---|
253 | and the solving of the elliptic equation associated with the surface pressure |
---|
254 | gradient computation (delocalization over the whole horizontal domain). |
---|
255 | Therefore, a pencil strategy is used for the data sub-structuration |
---|
256 | \gmcomment{no idea what this means!} |
---|
257 | : the 3D initial domain is laid out on local processor |
---|
258 | memories following a 2D horizontal topological splitting. Each sub-domain |
---|
259 | computes its own surface and bottom boundary conditions and has a side |
---|
260 | wall overlapping interface which defines the lateral boundary conditions for |
---|
261 | computations in the inner sub-domain. The overlapping area consists of the |
---|
262 | two rows at each edge of the sub-domain. After a computation, a communication |
---|
263 | phase starts: each processor sends to its neighbouring processors the update |
---|
264 | values of the points corresponding to the interior overlapping area to its |
---|
265 | neighbouring sub-domain (i.e. the innermost of the two overlapping rows). |
---|
266 | The communication is done through message passing. Usually the parallel virtual |
---|
267 | language, PVM, is used as it is a standard language available on nearly all |
---|
268 | MPP computers. More specific languages (i.e. computer dependant languages) |
---|
269 | can be easily used to speed up the communication, such as SHEM on a T3E |
---|
270 | computer. The data exchanges between processors are required at the very |
---|
271 | place where lateral domain boundary conditions are set in the mono-domain |
---|
272 | computation (\S III.10-c): the lbc\_lnk routine which manages such conditions |
---|
273 | is substituted by mpplnk.F or mpplnk2.F routine when running on an MPP |
---|
274 | computer (\key{mpp\_mpi} defined). It has to be pointed out that when using |
---|
275 | the MPP version of the model, the east-west cyclic boundary condition is done |
---|
276 | implicitly, whilst the south-symmetric boundary condition option is not available. |
---|
277 | |
---|
278 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
279 | \begin{figure}[!t] \begin{center} |
---|
280 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mpp.pdf} |
---|
281 | \caption{ \label{Fig_mpp} |
---|
282 | Positioning of a sub-domain when massively parallel processing is used. } |
---|
283 | \end{center} \end{figure} |
---|
284 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
285 | |
---|
286 | In the standard version of the OPA model, the splitting is regular and arithmetic. |
---|
287 | the i-axis is divided by \jp{jpni} and the j-axis by \jp{jpnj} for a number of processors |
---|
288 | \jp{jpnij} most often equal to $jpni \times jpnj$ (model parameters set in |
---|
289 | \mdl{par\_oce}). Each processor is independent and without message passing |
---|
290 | or synchronous process |
---|
291 | \gmcomment{how does a synchronous process relate to this?}, |
---|
292 | programs run alone and access just its own local memory. For this reason, the |
---|
293 | main model dimensions are now the local dimensions of the subdomain (pencil) |
---|
294 | that are named \jp{jpi}, \jp{jpj}, \jp{jpk}. These dimensions include the internal |
---|
295 | domain and the overlapping rows. The number of rows to exchange (known as |
---|
296 | the halo) is usually set to one (\jp{jpreci}=1, in \mdl{par\_oce}). The whole domain |
---|
297 | dimensions are named \jp{jpiglo}, \jp{jpjglo} and \jp{jpk}. The relationship between |
---|
298 | the whole domain and a sub-domain is: |
---|
299 | \begin{eqnarray} |
---|
300 | jpi & = & ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci \nonumber \\ |
---|
301 | jpj & = & ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj \label{Eq_lbc_jpi} |
---|
302 | \end{eqnarray} |
---|
303 | where \jp{jpni}, \jp{jpnj} are the number of processors following the i- and j-axis. |
---|
304 | |
---|
305 | \colorbox{yellow}{Figure IV.3: example of a domain splitting with 9 processors and |
---|
306 | no east-west cyclic boundary conditions.} |
---|
307 | |
---|
308 | One also defines variables nldi and nlei which correspond to the internal |
---|
309 | domain bounds, and the variables nimpp and njmpp which are the position |
---|
310 | of the (1,1) grid-point in the global domain. An element of $T_{l}$, a local array |
---|
311 | (subdomain) corresponds to an element of $T_{g}$, a global array |
---|
312 | (whole domain) by the relationship: |
---|
313 | \begin{equation} \label{Eq_lbc_nimpp} |
---|
314 | T_{g} (i+nimpp-1,j+njmpp-1,k) = T_{l} (i,j,k), |
---|
315 | \end{equation} |
---|
316 | with $1 \leq i \leq jpi$, $1 \leq j \leq jpj $ , and $1 \leq k \leq jpk$. |
---|
317 | |
---|
318 | Processors are numbered from 0 to $jpnij-1$, the number is saved in the variable |
---|
319 | nproc. In the standard version, a processor has no more than four neighbouring |
---|
320 | processors named nono (for north), noea (east), noso (south) and nowe (west) |
---|
321 | and two variables, nbondi and nbondj, indicate the relative position of the processor |
---|
322 | \colorbox{yellow}{(see Fig.IV.3)}: |
---|
323 | \begin{itemize} |
---|
324 | \item nbondi = -1 an east neighbour, no west processor, |
---|
325 | \item nbondi = 0 an east neighbour, a west neighbour, |
---|
326 | \item nbondi = 1 no east processor, a west neighbour, |
---|
327 | \item nbondi = 2 no splitting following the i-axis. |
---|
328 | \end{itemize} |
---|
329 | During the simulation, processors exchange data with their neighbours. |
---|
330 | If there is effectively a neighbour, the processor receives variables from this |
---|
331 | processor on its overlapping row, and sends the data issued from internal |
---|
332 | domain corresponding to the overlapping row of the other processor. |
---|
333 | |
---|
334 | \colorbox{yellow}{Figure IV.4: pencil splitting with the additional outer halos } |
---|
335 | |
---|
336 | |
---|
337 | The \NEMO model computes equation terms with the help of mask arrays (0 on land |
---|
338 | points and 1 on sea points). It is easily readable and very efficient in the context of |
---|
339 | a computer with vectorial architecture. However, in the case of a scalar processor, |
---|
340 | computations over the land regions become more expensive in terms of CPU time. |
---|
341 | It is worse when we use a complex configuration with a realistic bathymetry like the |
---|
342 | global ocean where more than 50 \% of points are land points. For this reason, a |
---|
343 | pre-processing tool can be used to choose the mpp domain decomposition with a |
---|
344 | maximum number of only land points processors, which can then be eliminated. |
---|
345 | (For example, the mpp\_optimiz tools, available from the DRAKKAR web site.) |
---|
346 | This optimisation is dependent on the specific bathymetry employed. The user |
---|
347 | then chooses optimal parameters \jp{jpni}, \jp{jpnj} and \jp{jpnij} with |
---|
348 | $jpnij < jpni \times jpnj$, leading to the elimination of $jpni \times jpnj - jpnij$ |
---|
349 | land processors. When those parameters are specified in module \mdl{par\_oce}, |
---|
350 | the algorithm in the \rou{inimpp2} routine sets each processor's parameters (nbound, |
---|
351 | nono, noea,...) so that the land-only processors are not taken into account. |
---|
352 | |
---|
353 | \colorbox{yellow}{Note that the inimpp2 routine is general so that the original inimpp |
---|
354 | routine should be suppressed from the code.} |
---|
355 | |
---|
356 | When land processors are eliminated, the value corresponding to these locations in |
---|
357 | the model output files is zero. Note that this is a problem for a mesh output file written |
---|
358 | by such a model configuration, because model users often divide by the scale factors |
---|
359 | ($e1t$, $e2t$, etc) and do not expect the grid size to be zero, even on land. It may be |
---|
360 | best not to eliminate land processors when running the model especially to write the |
---|
361 | mesh files as outputs (when \np{nn\_msh} namelist parameter differs from 0). |
---|
362 | %% |
---|
363 | \gmcomment{Steven : dont understand this, no land processor means no output file |
---|
364 | covering this part of globe; its only when files are stitched together into one that you |
---|
365 | can leave a hole} |
---|
366 | %% |
---|
367 | |
---|
368 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
369 | \begin{figure}[!ht] \begin{center} |
---|
370 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_mppini2.pdf} |
---|
371 | \caption { \label{Fig_mppini2} |
---|
372 | Example of Atlantic domain defined for the CLIPPER projet. Initial grid is |
---|
373 | composed of 773 x 1236 horizontal points. |
---|
374 | (a) the domain is split onto 9 \time 20 subdomains (jpni=9, jpnj=20). |
---|
375 | 52 subdomains are land areas. |
---|
376 | (b) 52 subdomains are eliminated (white rectangles) and the resulting number |
---|
377 | of processors really used during the computation is jpnij=128.} |
---|
378 | \end{center} \end{figure} |
---|
379 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
380 | |
---|
381 | |
---|
382 | % ================================================================ |
---|
383 | % Open Boundary Conditions |
---|
384 | % ================================================================ |
---|
385 | \section{Open Boundary Conditions (\key{obc}) (OBC)} |
---|
386 | \label{LBC_obc} |
---|
387 | %-----------------------------------------nam_obc ------------------------------------------- |
---|
388 | %- nobc_dta = 0 ! = 0 the obc data are equal to the initial state |
---|
389 | %- ! = 1 the obc data are read in 'obc .dta' files |
---|
390 | %- rn_dpein = 1. ! ??? |
---|
391 | %- rn_dpwin = 1. ! ??? |
---|
392 | %- rn_dpnin = 30. ! ??? |
---|
393 | %- rn_dpsin = 1. ! ??? |
---|
394 | %- rn_dpeob = 1500. ! time relaxation (days) for the east open boundary |
---|
395 | %- rn_dpwob = 15. ! " " for the west open boundary |
---|
396 | %- rn_dpnob = 150. ! " " for the north open boundary |
---|
397 | %- rn_dpsob = 15. ! " " for the south open boundary |
---|
398 | %- ln_obc_clim = .true. ! climatological obc data files (default T) |
---|
399 | %- ln_vol_cst = .true. ! total volume conserved |
---|
400 | \namdisplay{namobc} |
---|
401 | |
---|
402 | It is often necessary to implement a model configuration limited to an oceanic |
---|
403 | region or a basin, which communicates with the rest of the global ocean through |
---|
404 | ''open boundaries''. As stated by \citet{Roed1986}, an open boundary is a |
---|
405 | computational border where the aim of the calculations is to allow the perturbations |
---|
406 | generated inside the computational domain to leave it without deterioration of the |
---|
407 | inner model solution. However, an open boundary also has to let information from |
---|
408 | the outer ocean enter the model and should support inflow and outflow conditions. |
---|
409 | |
---|
410 | The open boundary package OBC is the first open boundary option developed in |
---|
411 | NEMO (originally in OPA8.2). It allows the user to |
---|
412 | \begin{itemize} |
---|
413 | \item tell the model that a boundary is ''open'' and not closed by a wall, for example |
---|
414 | by modifying the calculation of the divergence of velocity there; |
---|
415 | \item impose values of tracers and velocities at that boundary (values which may |
---|
416 | be taken from a climatology): this is the``fixed OBC'' option. |
---|
417 | \item calculate boundary values by a sophisticated algorithm combining radiation |
---|
418 | and relaxation (``radiative OBC'' option) |
---|
419 | \end{itemize} |
---|
420 | |
---|
421 | The package resides in the OBC directory. It is described here in four parts: the |
---|
422 | boundary geometry (parameters to be set in \mdl{obc\_par}), the forcing data at |
---|
423 | the boundaries (module \mdl{obcdta}), the radiation algorithm involving the |
---|
424 | namelist and module \mdl{obcrad}, and a brief presentation of boundary update |
---|
425 | and restart files. |
---|
426 | |
---|
427 | %---------------------------------------------- |
---|
428 | \subsection{Boundary geometry} |
---|
429 | \label{OBC_geom} |
---|
430 | % |
---|
431 | First one has to realize that open boundaries may not necessarily be located |
---|
432 | at the extremities of the computational domain. They may exist in the middle |
---|
433 | of the domain, for example at Gibraltar Straits if one wants to avoid including |
---|
434 | the Mediterranean in an Atlantic domain. This flexibility has been found necessary |
---|
435 | for the CLIPPER project \citep{Treguier_al_JGR01}. Because of the complexity of the |
---|
436 | geometry of ocean basins, it may even be necessary to have more than one |
---|
437 | ''west'' open boundary, more than one ''north'', etc. This is not possible with |
---|
438 | the OBC option: only one open boundary of each kind, west, east, south and |
---|
439 | north is allowed; these names refer to the grid geometry (not to the direction |
---|
440 | of the geographical ''west'', ''east'', etc). |
---|
441 | |
---|
442 | The open boundary geometry is set by a series of parameters in the module |
---|
443 | \mdl{obc\_par}. For an eastern open boundary, parameters are \jp{lp\_obc\_east} |
---|
444 | (true if an east open boundary exists), \jp{jpieob} the $i$-index along which |
---|
445 | the eastern open boundary (eob) is located, \jp{jpjed} the $j$-index at which |
---|
446 | it starts, and \jp{jpjef} the $j$-index where it ends (note $d$ is for ''d\'{e}but'' |
---|
447 | and $f$ for ''fin'' in French). Similar parameters exist for the west, south and |
---|
448 | north cases (Table~\ref{Tab_obc_param}). |
---|
449 | |
---|
450 | |
---|
451 | %--------------------------------------------------TABLE-------------------------------------------------- |
---|
452 | \begin{table}[htbp] \begin{center} \begin{tabular}{|l|c|c|c|} |
---|
453 | \hline |
---|
454 | Boundary and & Constant index & Starting index (d\'{e}but) & Ending index (fin) \\ |
---|
455 | Logical flag & & & \\ |
---|
456 | \hline |
---|
457 | West & \jp{jpiwob} $>= 2$ & \jp{jpjwd}$>= 2$ & \jp{jpjwf}<= \jp{jpjglo}-1 \\ |
---|
458 | lp\_obc\_west & $i$-index of a $u$ point & $j$ of a $T$ point &$j$ of a $T$ point \\ |
---|
459 | \hline |
---|
460 | East & \jp{jpieob}$<=$\jp{jpiglo}-2&\jp{jpjed} $>= 2$ & \jp{jpjef}$<=$ \jp{jpjglo}-1 \\ |
---|
461 | lp\_obc\_east & $i$-index of a $u$ point & $j$ of a $T$ point & $j$ of a $T$ point \\ |
---|
462 | \hline |
---|
463 | South & \jp{jpjsob} $>= 2$ & \jp{jpisd} $>= 2$ & \jp{jpisf}$<=$\jp{jpiglo}-1 \\ |
---|
464 | lp\_obc\_south & $j$-index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\ |
---|
465 | \hline |
---|
466 | North & \jp{jpjnob} $<=$ \jp{jpjglo}-2& \jp{jpind} $>= 2$ & \jp{jpinf}$<=$\jp{jpiglo}-1 \\ |
---|
467 | lp\_obc\_north & $j$-index of a $v$ point & $i$ of a $T$ point & $i$ of a $T$ point \\ |
---|
468 | \hline |
---|
469 | \end{tabular} \end{center} |
---|
470 | \caption{ \label{Tab_obc_param} |
---|
471 | Names of different indices relating to the open boundaries. In the case |
---|
472 | of a completely open ocean domain with four ocean boundaries, the parameters |
---|
473 | take exactly the values indicated.} |
---|
474 | \end{table} |
---|
475 | %------------------------------------------------------------------------------------------------------------ |
---|
476 | |
---|
477 | The open boundaries must be along coordinate lines. On the C-grid, the boundary |
---|
478 | itself is along a line of normal velocity points: $v$ points for a zonal open boundary |
---|
479 | (the south or north one), and $u$ points for a meridional open boundary (the west |
---|
480 | or east one). Another constraint is that there still must be a row of masked points |
---|
481 | all around the domain, as if the domain were a closed basin (unless periodic conditions |
---|
482 | are used together with open boundary conditions). Therefore, an open boundary |
---|
483 | cannot be located at the first/last index, namely, 1, \jp{jpiglo} or \jp{jpjglo}. Also, |
---|
484 | the open boundary algorithm involves calculating the normal velocity points situated |
---|
485 | just on the boundary, as well as the tangential velocity and temperature and salinity |
---|
486 | just outside the boundary. This means that for a west/south boundary, normal |
---|
487 | velocities and temperature are calculated at the same index \jp{jpiwob} and |
---|
488 | \jp{jpjsob}, respectively. For an east/north boundary, the normal velocity is |
---|
489 | calculated at index \jp{jpieob} and \jp{jpjnob}, but the ``outside'' temperature is |
---|
490 | at index \jp{jpieob}+1 and \jp{jpjnob}+1. This means that \jp{jpieob}, \jp{jpjnob} |
---|
491 | cannot be bigger than \jp{jpiglo}-2, \jp{jpjglo}-2. |
---|
492 | |
---|
493 | |
---|
494 | The starting and ending indices are to be thought of as $T$ point indices: in many |
---|
495 | cases they indicate the first land $T$-point, at the extremity of an open boundary |
---|
496 | (the coast line follows the $f$ grid points, see Fig.~\ref{Fig_obc_north} for an example |
---|
497 | of a northern open boundary). All indices are relative to the global domain. In the |
---|
498 | free surface case it is possible to have ``ocean corners'', that is, an open boundary |
---|
499 | starting and ending in the ocean. |
---|
500 | |
---|
501 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
502 | \begin{figure}[!t] \begin{center} |
---|
503 | \includegraphics[width=0.70\textwidth]{./TexFiles/Figures/Fig_obc_north.pdf} |
---|
504 | \caption{ \label{Fig_obc_north} |
---|
505 | Localization of the North open boundary points.} |
---|
506 | \end{center} \end{figure} |
---|
507 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
508 | |
---|
509 | Although not compulsory, it is highly recommended that the bathymetry in the |
---|
510 | vicinity of an open boundary follows the following rule: in the direction perpendicular |
---|
511 | to the open line, the water depth should be constant for 4 grid points. This is in |
---|
512 | order to ensure that the radiation condition, which involves model variables next |
---|
513 | to the boundary, is calculated in a consistent way. On Fig.\ref{Fig_obc_north} we |
---|
514 | indicate by an $=$ symbol, the points which should have the same depth. It means |
---|
515 | that at the 4 points near the boundary, the bathymetry is cylindrical \gmcomment{not sure |
---|
516 | why cylindrical}. The line behind the open $T$-line must be 0 in the bathymetry file |
---|
517 | (as shown on Fig.\ref{Fig_obc_north} for example). |
---|
518 | |
---|
519 | %---------------------------------------------- |
---|
520 | \subsection{Boundary data} |
---|
521 | \label{OBC_data} |
---|
522 | |
---|
523 | It is necessary to provide information at the boundaries. The simplest case is |
---|
524 | when this information does not change in time and is equal to the initial conditions |
---|
525 | (namelist variable \np{nn\_obcdta}=0). This is the case for the standard configuration |
---|
526 | EEL5 with open boundaries. When (\np{nn\_obcdta}=1), open boundary information |
---|
527 | is read from netcdf files. For convenience the input files are supposed to be similar |
---|
528 | to the ''history'' NEMO output files, for dimension names and variable names. |
---|
529 | Open boundary arrays must be dimensioned according to the parameters of table~ |
---|
530 | \ref{Tab_obc_param}: for example, at the western boundary, arrays have a |
---|
531 | dimension of \jp{jpwf}-\jp{jpwd}+1 in the horizontal and \jp{jpk} in the vertical. |
---|
532 | |
---|
533 | When ocean observations are used to generate the boundary data (a hydrographic |
---|
534 | section for example, as in \citet{Treguier_al_JGR01}) it happens often that only the velocity |
---|
535 | normal to the boundary is known, which is the reason why the initial OBC code |
---|
536 | assumes that only $T$, $S$, and the normal velocity ($u$ or $v$) needs to be |
---|
537 | specified. As more and more global model solutions and ocean analysis products |
---|
538 | become available, it will be possible to provide information about all the variables |
---|
539 | (including the tangential velocity) so that the specification of four variables at each |
---|
540 | boundaries will become standard. For the sea surface height, one must distinguish |
---|
541 | between the filtered free surface case and the time-splitting or explicit treatment of |
---|
542 | the free surface. |
---|
543 | In the first case, it is assumed that the user does not wish to represent high |
---|
544 | frequency motions such as tides. The boundary condition is thus one of zero |
---|
545 | normal gradient of sea surface height at the open boundaries, following \citet{Marchesiello2001}. |
---|
546 | No information other than the total velocity needs to be provided at the open |
---|
547 | boundaries in that case. In the other two cases (time splitting or explicit free surface), |
---|
548 | the user must provide barotropic information (sea surface height and barotropic |
---|
549 | velocities) and the use of the Flather algorithm for barotropic variables is |
---|
550 | recommanded. However, this algorithm has not yet been fully tested and bugs |
---|
551 | remain in NEMO v2.3. Users should read the code carefully before using it. Finally, |
---|
552 | in the case of the rigid lid approximation the barotropic streamfunction must be |
---|
553 | provided, as documented in \citet{Treguier_al_JGR01}). This option is no longer |
---|
554 | recommended but remains in NEMO V2.3. |
---|
555 | |
---|
556 | One frequently encountered case is when an open boundary domain is constructed |
---|
557 | from a global or larger scale NEMO configuration. Assuming the domain corresponds |
---|
558 | to indices $ib:ie$, $jb:je$ of the global domain, the bathymetry and forcing of the |
---|
559 | small domain can be created by using the following netcdf utility on the global files: |
---|
560 | ncks -F $-d\;x,ib,ie$ $-d\;y,jb,je$ (part of the nco series of utilities, |
---|
561 | see their \href{http://nco.sourceforge.net}{website}). |
---|
562 | The open boundary files can be constructed using ncks |
---|
563 | commands, following table~\ref{Tab_obc_ind}. |
---|
564 | |
---|
565 | %--------------------------------------------------TABLE-------------------------------------------------- |
---|
566 | \begin{table}[htbp] \begin{center} \begin{tabular}{|l|c|c|c|c|c|} |
---|
567 | \hline |
---|
568 | OBC & Variable & file name & Index & Start & end \\ |
---|
569 | West & T,S & obcwest\_TS.nc & $ib$+1 & $jb$+1 & $je-1$ \\ |
---|
570 | & U & obcwest\_U.nc & $ib$+1 & $jb$+1 & $je-1$ \\ |
---|
571 | & V & obcwest\_V.nc & $ib$+1 & $jb$+1 & $je-1$ \\ |
---|
572 | \hline |
---|
573 | East & T,S & obceast\_TS.nc & $ie$-1 & $jb$+1 & $je-1$ \\ |
---|
574 | & U & obceast\_U.nc & $ie$-2 & $jb$+1 & $je-1$ \\ |
---|
575 | & V & obceast\_V.nc & $ie$-1 & $jb$+1 & $je-1$ \\ |
---|
576 | \hline |
---|
577 | South & T,S & obcsouth\_TS.nc & $jb$+1 & $ib$+1 & $ie-1$ \\ |
---|
578 | & U & obcsouth\_U.nc & $jb$+1 & $ib$+1 & $ie-1$ \\ |
---|
579 | & V & obcsouth\_V.nc & $jb$+1 & $ib$+1 & $ie-1$ \\ |
---|
580 | \hline |
---|
581 | North & T,S & obcnorth\_TS.nc & $je$-1 & $ib$+1 & $ie-1$ \\ |
---|
582 | & U & obcnorth\_U.nc & $je$-1 & $ib$+1 & $ie-1$ \\ |
---|
583 | & V & obcnorth\_V.nc & $je$-2 & $ib$+1 & $ie-1$ \\ |
---|
584 | \hline |
---|
585 | \end{tabular} \end{center} |
---|
586 | \caption{ \label{Tab_obc_ind} |
---|
587 | Requirements for creating open boundary files from a global configuration, |
---|
588 | appropriate for the subdomain of indices $ib:ie$, $jb:je$. ``Index'' designates the |
---|
589 | $i$ or $j$ index along which the $u$ of $v$ boundary point is situated in the global |
---|
590 | configuration, starting and ending with the $j$ or $i$ indices indicated. |
---|
591 | For example, to generate file obcnorth\_V.nc, use the command ncks |
---|
592 | $-F$ $-d\;y,je-2$ $-d\;x,ib+1,ie-1$ } |
---|
593 | \end{table} |
---|
594 | %----------------------------------------------------------------------------------------------------------- |
---|
595 | |
---|
596 | It is assumed that the open boundary files contain the variables for the period of |
---|
597 | the model integration. If the boundary files contain one time frame, the boundary |
---|
598 | data is held fixed in time. If the files contain 12 values, it is assumed that the input |
---|
599 | is a climatology for a repeated annual cycle (corresponding to the case \np{ln\_obc\_clim} |
---|
600 | =true). The case of an arbitrary number of time frames is not yet implemented |
---|
601 | correctly; the user is required to write his own code in the module \mdl{obc\_dta} |
---|
602 | to deal with this situation. |
---|
603 | |
---|
604 | \subsection{Radiation algorithm} |
---|
605 | \label{OBC_rad} |
---|
606 | |
---|
607 | The art of open boundary management consists in applying a constraint strong |
---|
608 | enough that the inner domain "feels" the rest of the ocean, but weak enough |
---|
609 | that perturbations are allowed to leave the domain with minimum false reflections |
---|
610 | of energy. The constraints are specified separately at each boundary as time |
---|
611 | scales for ''inflow'' and ''outflow'' as defined below. The time scales are set (in days) |
---|
612 | by namelist parameters such as \np{rn\_dpein}, \np{rn\_dpeob} for the eastern open |
---|
613 | boundary for example. When both time scales are zero for a given boundary |
---|
614 | ($e.g.$ for the western boundary, \jp{lp\_obc\_west}=true, \np{rn\_dpwob}=0 and |
---|
615 | \np{rn\_dpwin}=0) this means that the boundary in question is a ''fixed '' boundary |
---|
616 | where the solution is set exactly by the boundary data. This is not recommended, |
---|
617 | except in combination with increased viscosity in a ''sponge'' layer next to the |
---|
618 | boundary in order to avoid spurious reflections. |
---|
619 | |
---|
620 | |
---|
621 | The radiation\/relaxation \gmcomment{the / doesnt seem to appear in the output} |
---|
622 | algorithm is applied when either relaxation time (for ''inflow'' or ''outflow'') is |
---|
623 | non-zero. It has been developed and tested in the SPEM model and its |
---|
624 | successor ROMS \citep{Barnier1996, Marchesiello2001}, which is an |
---|
625 | $s$-coordinate model on an Arakawa C-grid. Although the algorithm has |
---|
626 | been numerically successful in the CLIPPER Atlantic models, the physics |
---|
627 | do not work as expected \citep{Treguier_al_JGR01}. Users are invited to consider |
---|
628 | open boundary conditions (OBC hereafter) with some scepticism |
---|
629 | \citep{Durran2001, Blayo2005}. |
---|
630 | |
---|
631 | The first part of the algorithm calculates a phase velocity to determine |
---|
632 | whether perturbations tend to propagate toward, or away from, the |
---|
633 | boundary. Let us consider a model variable $\phi$. |
---|
634 | The phase velocities ($C_{\phi x}$,$C_{\phi y}$) for the variable $\phi$, |
---|
635 | in the directions normal and tangential to the boundary are |
---|
636 | \begin{equation} \label{Eq_obc_cphi} |
---|
637 | C_{\phi x} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{x} |
---|
638 | \;\;\;\;\; \;\;\; |
---|
639 | C_{\phi y} = \frac{ -\phi_{t} }{ ( \phi_{x}^{2} + \phi_{y}^{2}) } \phi_{y}. |
---|
640 | \end{equation} |
---|
641 | Following \citet{Treguier_al_JGR01} and \citet{Marchesiello2001} we retain only |
---|
642 | the normal component of the velocity, $C_{\phi x}$, setting $C_{\phi y} =0$ |
---|
643 | (but unlike the original Orlanski radiation algorithm we retain $\phi_{y}$ in |
---|
644 | the expression for $C_{\phi x}$). |
---|
645 | |
---|
646 | The discrete form of (\ref{Eq_obc_cphi}), described by \citet{Barnier1998}, |
---|
647 | takes into account the two rows of grid points situated inside the domain |
---|
648 | next to the boundary, and the three previous time steps ($n$, $n-1$, |
---|
649 | and $n-2$). The same equation can then be discretized at the boundary at |
---|
650 | time steps $n-1$, $n$ and $n+1$ \gmcomment{since the original was three time-level} |
---|
651 | in order to extrapolate for the new boundary value $\phi^{n+1}$. |
---|
652 | |
---|
653 | In the open boundary algorithm as implemented in NEMO v2.3, the new boundary |
---|
654 | values are updated differently depending on the sign of $C_{\phi x}$. Let us take |
---|
655 | an eastern boundary as an example. The solution for variable $\phi$ at the |
---|
656 | boundary is given by a generalized wave equation with phase velocity $C_{\phi}$, |
---|
657 | with the addition of a relaxation term, as: |
---|
658 | \begin{eqnarray} |
---|
659 | \phi_{t} & = & -C_{\phi x} \phi_{x} + \frac{1}{\tau_{o}} (\phi_{c}-\phi) |
---|
660 | \;\;\; \;\;\; \;\;\; (C_{\phi x} > 0), \label{Eq_obc_rado} \\ |
---|
661 | \phi_{t} & = & \frac{1}{\tau_{i}} (\phi_{c}-\phi) |
---|
662 | \;\;\; \;\;\; \;\;\;\;\;\; (C_{\phi x} < 0), \label{Eq_obc_radi} |
---|
663 | \end{eqnarray} |
---|
664 | where $\phi_{c}$ is the estimate of $\phi$ at the boundary, provided as boundary |
---|
665 | data. Note that in (\ref{Eq_obc_rado}), $C_{\phi x}$ is bounded by the ratio |
---|
666 | $\delta x/\delta t$ for stability reasons. When $C_{\phi x}$ is eastward (outward |
---|
667 | propagation), the radiation condition (\ref{Eq_obc_rado}) is used. |
---|
668 | When $C_{\phi x}$ is westward (inward propagation), (\ref{Eq_obc_radi}) is |
---|
669 | used with a strong relaxation to climatology (usually $\tau_{i}=\np{rn\_dpein}=$1~day). |
---|
670 | Equation (\ref{Eq_obc_radi}) is solved with a Euler time-stepping scheme. As a |
---|
671 | consequence, setting $\tau_{i}$ smaller than, or equal to the time step is equivalent |
---|
672 | to a fixed boundary condition. A time scale of one day is usually a good compromise |
---|
673 | which guarantees that the inflow conditions remain close to climatology while ensuring |
---|
674 | numerical stability. |
---|
675 | |
---|
676 | In the case of a western boundary located in the Eastern Atlantic, \citet{Penduff_al_JGR00} |
---|
677 | have been able to implement the radiation algorithm without any boundary data, |
---|
678 | using persistence from the previous time step instead. This solution has not worked |
---|
679 | in other cases \citep{Treguier_al_JGR01}, so that the use of boundary data is recommended. |
---|
680 | Even in the outflow condition (\ref{Eq_obc_rado}), we have found it desirable to |
---|
681 | maintain a weak relaxation to climatology. The time step is usually chosen so as to |
---|
682 | be larger than typical turbulent scales (of order 1000~days \gmcomment{or maybe seconds?}). |
---|
683 | |
---|
684 | The radiation condition is applied to the model variables: temperature, salinity, |
---|
685 | tangential and normal velocities. For normal and tangential velocities, $u$ and $v$, |
---|
686 | radiation is applied with phase velocities calculated from $u$ and $v$ respectively. |
---|
687 | For the radiation of tracers, we use the phase velocity calculated from the tangential |
---|
688 | velocity in order to avoid calculating too many independent radiation velocities and |
---|
689 | because tangential velocities and tracers have the same position along the boundary |
---|
690 | on a C-grid. |
---|
691 | |
---|
692 | \subsection{Domain decomposition (\key{mpp\_mpi})} |
---|
693 | \label{OBC_mpp} |
---|
694 | When \key{mpp\_mpi} is active in the code, the computational domain is divided |
---|
695 | into rectangles that are attributed each to a different processor. The open boundary |
---|
696 | code is ``mpp-compatible'' up to a certain point. The radiation algorithm will not |
---|
697 | work if there is an mpp subdomain boundary parallel to the open boundary at the |
---|
698 | index of the boundary, or the grid point after (outside), or three grid points before |
---|
699 | (inside). On the other hand, there is no problem if an mpp subdomain boundary |
---|
700 | cuts the open boundary perpendicularly. These geometrical limitations must be |
---|
701 | checked for by the user (there is no safeguard in the code). |
---|
702 | The general principle for the open boundary mpp code is that loops over the open |
---|
703 | boundaries {not sure what this means} are performed on local indices (nie0, |
---|
704 | nie1, nje0, nje1 for an eastern boundary for instance) that are initialized in module |
---|
705 | \mdl{obc\_ini}. Those indices have relevant values on the processors that contain |
---|
706 | a segment of an open boundary. For processors that do not include an open |
---|
707 | boundary segment, the indices are such that the calculations within the loops are |
---|
708 | not performed. |
---|
709 | \gmcomment{I dont understand most of the last few sentences} |
---|
710 | |
---|
711 | Arrays of climatological data that are read from files are seen by all processors |
---|
712 | and have the same dimensions for all (for instance, for the eastern boundary, |
---|
713 | uedta(jpjglo,jpk,2)). On the other hand, the arrays for the calculation of radiation |
---|
714 | are local to each processor (uebnd(jpj,jpk,3,3) for instance). This allowed the |
---|
715 | CLIPPER model for example, to save on memory where the eastern boundary |
---|
716 | crossed 8 processors so that \jp{jpj} was much smaller than (\jp{jpjef}-\jp{jpjed}+1). |
---|
717 | |
---|
718 | \subsection{Volume conservation} |
---|
719 | \label{OBC_vol} |
---|
720 | |
---|
721 | It is necessary to control the volume inside a domain when using open boundaries. |
---|
722 | With fixed boundaries, it is enough to ensure that the total inflow/outflow has |
---|
723 | reasonable values (either zero or a value compatible with an observed volume |
---|
724 | balance). When using radiative boundary conditions it is necessary to have a |
---|
725 | volume constraint because each open boundary works independently from the |
---|
726 | others. The methodology used to control this volume is identical to the one |
---|
727 | coded in the ROMS model \citep{Marchesiello2001}. |
---|
728 | |
---|
729 | |
---|
730 | %---------------------------------------- EXTRAS |
---|
731 | \colorbox{yellow}{Explain obc\_vol{\ldots}} |
---|
732 | |
---|
733 | \colorbox{yellow}{OBC algorithm for update, OBC restart, list of routines where obc key appears{\ldots}} |
---|
734 | |
---|
735 | \colorbox{yellow}{OBC rigid lid? {\ldots}} |
---|
736 | |
---|
737 | % ==================================================================== |
---|
738 | % Unstructured open boundaries BDY |
---|
739 | % ==================================================================== |
---|
740 | \section{Unstructured Open Boundary Conditions (\key{bdy}) (BDY)} |
---|
741 | \label{LBC_bdy} |
---|
742 | |
---|
743 | %-----------------------------------------nambdy-------------------------------------------- |
---|
744 | \namdisplay{nambdy} |
---|
745 | %----------------------------------------------------------------------------------------------- |
---|
746 | %-----------------------------------------nambdy_index-------------------------------------------- |
---|
747 | \namdisplay{nambdy_index} |
---|
748 | %----------------------------------------------------------------------------------------------- |
---|
749 | %-----------------------------------------nambdy_dta-------------------------------------------- |
---|
750 | \namdisplay{nambdy_dta} |
---|
751 | %----------------------------------------------------------------------------------------------- |
---|
752 | %-----------------------------------------nambdy_dta-------------------------------------------- |
---|
753 | \namdisplay{nambdy_dta2} |
---|
754 | %----------------------------------------------------------------------------------------------- |
---|
755 | |
---|
756 | The BDY module is an alternative implementation of open boundary |
---|
757 | conditions for regional configurations. It implements the Flow |
---|
758 | Relaxation Scheme algorithm for temperature, salinity, velocities and |
---|
759 | ice fields, and the Flather radiation condition for the depth-mean |
---|
760 | transports. The specification of the location of the open boundary is |
---|
761 | completely flexible and allows for example the open boundary to follow |
---|
762 | an isobath or other irregular contour. |
---|
763 | |
---|
764 | The BDY module was modelled on the OBC module and shares many features |
---|
765 | and a similar coding structure \citep{Chanut2005}. |
---|
766 | |
---|
767 | The BDY module is completely rewritten at NEMO 3.4 and there is a new |
---|
768 | set of namelists. Boundary data files used with earlier versions of |
---|
769 | NEMO may need to be re-ordered to work with this version. See the |
---|
770 | section on the Input Boundary Data Files for details. |
---|
771 | |
---|
772 | %---------------------------------------------- |
---|
773 | \subsection{The namelists} |
---|
774 | \label{BDY_namelist} |
---|
775 | |
---|
776 | It is possible to define more than one boundary ``set'' and apply |
---|
777 | different boundary conditions to each set. The number of boundary |
---|
778 | sets is defined by \np{nb\_bdy}. Each boundary set may be defined |
---|
779 | as a set of straight line segments in a namelist |
---|
780 | (\np{ln\_coords\_file}=.false.) or read in from a file |
---|
781 | (\np{ln\_coords\_file}=.true.). If the set is defined in a namelist, |
---|
782 | then the namelists nambdy\_index must be included separately, one for |
---|
783 | each set. If the set is defined by a file, then a |
---|
784 | ``coordinates.bdy.nc'' file must be provided. The coordinates.bdy file |
---|
785 | is analagous to the usual NEMO ``coordinates.nc'' file. In the example |
---|
786 | above, there are two boundary sets, the first of which is defined via |
---|
787 | a file and the second is defined in a namelist. For more details of |
---|
788 | the definition of the boundary geometry see section |
---|
789 | \ref{BDY_geometry}. |
---|
790 | |
---|
791 | For each boundary set a boundary |
---|
792 | condition has to be chosen for the barotropic solution (``u2d'': |
---|
793 | sea-surface height and barotropic velocities), for the baroclinic |
---|
794 | velocities (``u3d''), and for the active tracers\footnote{The BDY |
---|
795 | module does not deal with passive tracers at this version} |
---|
796 | (``tra''). For each set of variables there is a choice of algorithm |
---|
797 | and a choice for the data, eg. for the active tracers the algorithm is |
---|
798 | set by \np{nn\_tra} and the choice of data is set by |
---|
799 | \np{nn\_tra\_dta}. |
---|
800 | |
---|
801 | The choice of algorithm is currently as follows: |
---|
802 | |
---|
803 | \mbox{} |
---|
804 | |
---|
805 | \begin{itemize} |
---|
806 | \item[0.] No boundary condition applied. So the solution will ``see'' |
---|
807 | the land points around the edge of the edge of the domain. |
---|
808 | \item[1.] Flow Relaxation Scheme (FRS) available for all variables. |
---|
809 | \item[2.] Flather radiation scheme for the barotropic variables. The |
---|
810 | Flather scheme is not compatible with the filtered free surface |
---|
811 | ({\it dynspg\_ts}). |
---|
812 | \end{itemize} |
---|
813 | |
---|
814 | \mbox{} |
---|
815 | |
---|
816 | The main choice for the boundary data is |
---|
817 | to use initial conditions as boundary data (\np{nn\_tra\_dta}=0) or to |
---|
818 | use external data from a file (\np{nn\_tra\_dta}=1). For the |
---|
819 | barotropic solution there is also the option to use tidal |
---|
820 | harmonic forcing either by itself or in addition to other external |
---|
821 | data. |
---|
822 | |
---|
823 | If external boundary data is required then the nambdy\_dta namelist |
---|
824 | must be defined. One nambdy\_dta namelist is required for each boundary |
---|
825 | set in the order in which the boundary sets are defined in nambdy. In |
---|
826 | the example given, two boundary sets have been defined and so there |
---|
827 | are two nambdy\_dta namelists. The boundary data is read in using the |
---|
828 | fldread module, so the nambdy\_dta namelist is in the format required |
---|
829 | for fldread. For each variable required, the filename, the frequency |
---|
830 | of the files and the frequency of the data in the files is given. Also |
---|
831 | whether or not time-interpolation is required and whether the data is |
---|
832 | climatological (time-cyclic) data. Note that on-the-fly spatial |
---|
833 | interpolation of boundary data is not available at this version. |
---|
834 | |
---|
835 | In the example namelists given, two boundary sets are defined. The |
---|
836 | first set is defined via a file and applies FRS conditions to |
---|
837 | temperature and salinity and Flather conditions to the barotropic |
---|
838 | variables. External data is provided in daily files (from a |
---|
839 | large-scale model). Tidal harmonic forcing is also used. The second |
---|
840 | set is defined in a namelist. FRS conditions are applied on |
---|
841 | temperature and salinity and climatological data is read from external |
---|
842 | files. |
---|
843 | |
---|
844 | %---------------------------------------------- |
---|
845 | \subsection{The Flow Relaxation Scheme} |
---|
846 | \label{BDY_FRS_scheme} |
---|
847 | |
---|
848 | The Flow Relaxation Scheme (FRS) \citep{Davies_QJRMS76,Engerdahl_Tel95}, |
---|
849 | applies a simple relaxation of the model fields to |
---|
850 | externally-specified values over a zone next to the edge of the model |
---|
851 | domain. Given a model prognostic variable $\Phi$ |
---|
852 | \begin{equation} \label{Eq_bdy_frs1} |
---|
853 | \Phi(d) = \alpha(d)\Phi_{e}(d) + (1-\alpha(d))\Phi_{m}(d)\;\;\;\;\; d=1,N |
---|
854 | \end{equation} |
---|
855 | where $\Phi_{m}$ is the model solution and $\Phi_{e}$ is the specified |
---|
856 | external field, $d$ gives the discrete distance from the model |
---|
857 | boundary and $\alpha$ is a parameter that varies from $1$ at $d=1$ to |
---|
858 | a small value at $d=N$. It can be shown that this scheme is equivalent |
---|
859 | to adding a relaxation term to the prognostic equation for $\Phi$ of |
---|
860 | the form: |
---|
861 | \begin{equation} \label{Eq_bdy_frs2} |
---|
862 | -\frac{1}{\tau}\left(\Phi - \Phi_{e}\right) |
---|
863 | \end{equation} |
---|
864 | where the relaxation time scale $\tau$ is given by a function of |
---|
865 | $\alpha$ and the model time step $\Delta t$: |
---|
866 | \begin{equation} \label{Eq_bdy_frs3} |
---|
867 | \tau = \frac{1-\alpha}{\alpha} \,\rdt |
---|
868 | \end{equation} |
---|
869 | Thus the model solution is completely prescribed by the external |
---|
870 | conditions at the edge of the model domain and is relaxed towards the |
---|
871 | external conditions over the rest of the FRS zone. The application of |
---|
872 | a relaxation zone helps to prevent spurious reflection of outgoing |
---|
873 | signals from the model boundary. |
---|
874 | |
---|
875 | The function $\alpha$ is specified as a $tanh$ function: |
---|
876 | \begin{equation} \label{Eq_bdy_frs4} |
---|
877 | \alpha(d) = 1 - \tanh\left(\frac{d-1}{2}\right), \quad d=1,N |
---|
878 | \end{equation} |
---|
879 | The width of the FRS zone is specified in the namelist as |
---|
880 | \np{nn\_rimwidth}. This is typically set to a value between 8 and 10. |
---|
881 | |
---|
882 | %---------------------------------------------- |
---|
883 | \subsection{The Flather radiation scheme} |
---|
884 | \label{BDY_flather_scheme} |
---|
885 | |
---|
886 | The \citet{Flather_JPO94} scheme is a radiation condition on the normal, depth-mean |
---|
887 | transport across the open boundary. It takes the form |
---|
888 | \begin{equation} \label{Eq_bdy_fla1} |
---|
889 | U = U_{e} + \frac{c}{h}\left(\eta - \eta_{e}\right), |
---|
890 | \end{equation} |
---|
891 | where $U$ is the depth-mean velocity normal to the boundary and $\eta$ |
---|
892 | is the sea surface height, both from the model. The subscript $e$ |
---|
893 | indicates the same fields from external sources. The speed of external |
---|
894 | gravity waves is given by $c = \sqrt{gh}$, and $h$ is the depth of the |
---|
895 | water column. The depth-mean normal velocity along the edge of the |
---|
896 | model domain is set equal to the |
---|
897 | external depth-mean normal velocity, plus a correction term that |
---|
898 | allows gravity waves generated internally to exit the model boundary. |
---|
899 | Note that the sea-surface height gradient in \eqref{Eq_bdy_fla1} |
---|
900 | is a spatial gradient across the model boundary, so that $\eta_{e}$ is |
---|
901 | defined on the $T$ points with $nbr=1$ and $\eta$ is defined on the |
---|
902 | $T$ points with $nbr=2$. $U$ and $U_{e}$ are defined on the $U$ or |
---|
903 | $V$ points with $nbr=1$, $i.e.$ between the two $T$ grid points. |
---|
904 | |
---|
905 | %---------------------------------------------- |
---|
906 | \subsection{Boundary geometry} |
---|
907 | \label{BDY_geometry} |
---|
908 | |
---|
909 | Each open boundary set is defined as a list of points. The information |
---|
910 | is stored in the arrays $nbi$, $nbj$, and $nbr$ in the $idx\_bdy$ |
---|
911 | structure. The $nbi$ and $nbj$ arrays |
---|
912 | define the local $(i,j)$ indices of each point in the boundary zone |
---|
913 | and the $nbr$ array defines the discrete distance from the boundary |
---|
914 | with $nbr=1$ meaning that the point is next to the edge of the |
---|
915 | model domain and $nbr>1$ showing that the point is increasingly |
---|
916 | further away from the edge of the model domain. A set of $nbi$, $nbj$, |
---|
917 | and $nbr$ arrays is defined for each of the $T$, $U$ and $V$ |
---|
918 | grids. Figure \ref{Fig_LBC_bdy_geom} shows an example of an irregular |
---|
919 | boundary. |
---|
920 | |
---|
921 | The boundary geometry for each set may be defined in a namelist |
---|
922 | nambdy\_index or by reading in a ``coordinates.bdy.nc'' file. The |
---|
923 | nambdy\_index namelist defines a series of straight-line segments for |
---|
924 | north, east, south and west boundaries. For the northern boundary, |
---|
925 | \np{nbdysegn} gives the number of segments, \np{jpjnob} gives the $j$ |
---|
926 | index for each segment and \np{jpindt} and \np{jpinft} give the start |
---|
927 | and end $i$ indices for each segment with similar for the other |
---|
928 | boundaries. These segments define a list of $T$ grid points along the |
---|
929 | outermost row of the boundary ($nbr\,=\, 1$). The code deduces the $U$ and |
---|
930 | $V$ points and also the points for $nbr\,>\, 1$ if |
---|
931 | $nn\_rimwidth\,>\,1$. |
---|
932 | |
---|
933 | The boundary geometry may also be defined from a |
---|
934 | ``coordinates.bdy.nc'' file. Figure \ref{Fig_LBC_nc_header} |
---|
935 | gives an example of the header information from such a file. The file |
---|
936 | should contain the index arrays for each of the $T$, $U$ and $V$ |
---|
937 | grids. The arrays must be in order of increasing $nbr$. Note that the |
---|
938 | $nbi$, $nbj$ values in the file are global values and are converted to |
---|
939 | local values in the code. Typically this file will be used to generate |
---|
940 | external boundary data via interpolation and so will also contain the |
---|
941 | latitudes and longitudes of each point as shown. However, this is not |
---|
942 | necessary to run the model. |
---|
943 | |
---|
944 | For some choices of irregular boundary the model domain may contain |
---|
945 | areas of ocean which are not part of the computational domain. For |
---|
946 | example if an open boundary is defined along an isobath, say at the |
---|
947 | shelf break, then the areas of ocean outside of this boundary will |
---|
948 | need to be masked out. This can be done by reading a mask file defined |
---|
949 | as \np{cn\_mask\_file} in the nam\_bdy namelist. Only one mask file is |
---|
950 | used even if multiple boundary sets are defined. |
---|
951 | |
---|
952 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
953 | \begin{figure}[!t] \begin{center} |
---|
954 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_bdy_geom.pdf} |
---|
955 | \caption { \label{Fig_LBC_bdy_geom} |
---|
956 | Example of geometry of unstructured open boundary} |
---|
957 | \end{center} \end{figure} |
---|
958 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
959 | |
---|
960 | %---------------------------------------------- |
---|
961 | \subsection{Input boundary data files} |
---|
962 | \label{BDY_data} |
---|
963 | |
---|
964 | The data files contain the data arrays |
---|
965 | in the order in which the points are defined in the $nbi$ and $nbj$ |
---|
966 | arrays. The data arrays are dimensioned on: a time dimension; |
---|
967 | $xb$ which is the index of the boundary data point in the horizontal; |
---|
968 | and $yb$ which is a degenerate dimension of 1 to enable the file to be |
---|
969 | read by the standard NEMO I/O routines. The 3D fields also have a |
---|
970 | depth dimension. |
---|
971 | |
---|
972 | At Version 3.4 there are new restrictions on the order in which the |
---|
973 | boundary points are defined (and therefore restrictions on the order |
---|
974 | of the data in the file). In particular: |
---|
975 | |
---|
976 | \mbox{} |
---|
977 | |
---|
978 | \begin{enumerate} |
---|
979 | \item The data points must be in order of increasing $nbr$, ie. all |
---|
980 | the $nbr=1$ points, then all the $nbr=2$ points etc. |
---|
981 | \item All the data for a particular boundary set must be in the same |
---|
982 | order. (Prior to 3.4 it was possible to define barotropic data in a |
---|
983 | different order to the data for tracers and baroclinic velocities). |
---|
984 | \end{enumerate} |
---|
985 | |
---|
986 | \mbox{} |
---|
987 | |
---|
988 | These restrictions mean that data files used with previous versions of |
---|
989 | the model may not work with version 3.4. A fortran utility |
---|
990 | {\it bdy\_reorder} exists in the TOOLS directory which will re-order the |
---|
991 | data in old BDY data files. |
---|
992 | |
---|
993 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
994 | \begin{figure}[!t] \begin{center} |
---|
995 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_LBC_nc_header.pdf} |
---|
996 | \caption { \label{Fig_LBC_nc_header} |
---|
997 | Example of the header for a coordinates.bdy.nc file} |
---|
998 | \end{center} \end{figure} |
---|
999 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
1000 | |
---|
1001 | %---------------------------------------------- |
---|
1002 | \subsection{Volume correction} |
---|
1003 | \label{BDY_vol_corr} |
---|
1004 | |
---|
1005 | There is an option to force the total volume in the regional model to be constant, |
---|
1006 | similar to the option in the OBC module. This is controlled by the \np{nn\_volctl} |
---|
1007 | parameter in the namelist. A value of \np{nn\_volctl}~=~0 indicates that this option is not used. |
---|
1008 | If \np{nn\_volctl}~=~1 then a correction is applied to the normal velocities |
---|
1009 | around the boundary at each timestep to ensure that the integrated volume flow |
---|
1010 | through the boundary is zero. If \np{nn\_volctl}~=~2 then the calculation of |
---|
1011 | the volume change on the timestep includes the change due to the freshwater |
---|
1012 | flux across the surface and the correction velocity corrects for this as well. |
---|
1013 | |
---|
1014 | If more than one boundary set is used then volume correction is |
---|
1015 | applied to all boundaries at once. |
---|
1016 | |
---|
1017 | \newpage |
---|
1018 | %---------------------------------------------- |
---|
1019 | \subsection{Tidal harmonic forcing} |
---|
1020 | \label{BDY_tides} |
---|
1021 | |
---|
1022 | %-----------------------------------------nambdy_tide-------------------------------------------- |
---|
1023 | \namdisplay{nambdy_tide} |
---|
1024 | %----------------------------------------------------------------------------------------------- |
---|
1025 | |
---|
1026 | To be written.... |
---|
1027 | |
---|
1028 | |
---|
1029 | |
---|
1030 | |
---|