1 | \magnification =\magstep1 |
---|
2 | \count0=77 |
---|
3 | %definitions |
---|
4 | |
---|
5 | %end of definitions |
---|
6 | |
---|
7 | |
---|
8 | \centerline{ Olivier Thual, March 15$^{\rm th}$ 1992} |
---|
9 | \centerline{ Revised version June 30$^{\rm th}$ 1992} |
---|
10 | \bigskip |
---|
11 | |
---|
12 | \centerline{\bf SIMPLE OCEAN-ATMOSPHERE INTERPOLATION.} |
---|
13 | \centerline{\bf PART A: THE METHOD} |
---|
14 | |
---|
15 | \bigskip |
---|
16 | epi920629-naiva-ot6.tex |
---|
17 | A simple and naive method to interpolate fields between an atmospheric |
---|
18 | and an oceanic grid is presented. This crude method is coded into a set of |
---|
19 | FORTRAN programs, in the expectation of a more elaborated method which |
---|
20 | takes into account the particular physic of the ocean-atmosphere coupling. |
---|
21 | |
---|
22 | |
---|
23 | |
---|
24 | \beginsection 1. INTRODUCTION |
---|
25 | |
---|
26 | The coupling of an atmosphere General Circulation Model (GCM) and an |
---|
27 | oceanic GCM involves both a temporal coupling, which deals with |
---|
28 | the synchronization and the communication between the two codes at |
---|
29 | suitable time steps, and a spatial coupling which deals with the |
---|
30 | interpolations of suitable fields from one GCM grid to another, namely the |
---|
31 | fluxes from the atmospheric to the ocean grid, and the sea surface |
---|
32 | temperature (SST) in the reverse direction. |
---|
33 | |
---|
34 | |
---|
35 | |
---|
36 | These two aspects of the coupling are presently attacked in parallel in the |
---|
37 | CERFACS Climate Modelling \& Global Change team. They have started |
---|
38 | with a gathering of all the state of art informations that we have been |
---|
39 | able to record in the Climatic community and that have been summarized |
---|
40 | in [1]. Today, progress has been made in the two directions. |
---|
41 | |
---|
42 | In [2], a method based on the communication and synchronization of UNIX |
---|
43 | processes is developed for the GCMs temporal coupling task. The skeleton |
---|
44 | of a universal ``coupler'' has been written in FORTRAN. |
---|
45 | This software |
---|
46 | controls all the exchanges between an oceanic and an atmospheric codes. |
---|
47 | Its present conception allows a simple |
---|
48 | implementation in an already written and tested GCM, with minimum |
---|
49 | coding changes in it. |
---|
50 | |
---|
51 | |
---|
52 | |
---|
53 | This coupler skeleton includes a black box which will perform the spatial |
---|
54 | interpolations between the atmospheric and oceanic grids. Different |
---|
55 | interpolation schemes can be plugged in this black box, which further |
---|
56 | increases the universality of the coupler. The fundamental basis for the |
---|
57 | design of an advanced interpolation scheme is presented in [3]. In this |
---|
58 | scheme a particular focus is made on the ``local conservation'' of the |
---|
59 | fluxes when the interpolation from the atmospheric to the oceanic grid is |
---|
60 | performed. It is thought ([3]) than the non-conservation of the fluxes |
---|
61 | might be one of the causes of the drift problem encountered when two GCMs, |
---|
62 | individually well tuned, are coupled together. |
---|
63 | |
---|
64 | |
---|
65 | |
---|
66 | In the expectation of the achievement and coding of this advanced |
---|
67 | interpolation scheme, it has appeared necessary to plug into the coupler |
---|
68 | spatial interpolation black box a simple and quickly written package of |
---|
69 | FORTRAN subroutines in order to start real size tests of the coupler with |
---|
70 | GCMs. The purpose of this letter is to expose the |
---|
71 | ``naive'' and crude interpolation method that we have chosen in this spirit. |
---|
72 | It is based on the informations that we have collected in our inquiry on the |
---|
73 | state of the art spread among the climatic community (see [4]) and |
---|
74 | summarized in [1]. |
---|
75 | |
---|
76 | A basic introduction to the interpolation of data grid points is presented |
---|
77 | in Section 1, and the specific problem of the grid-to-grid interpolation is |
---|
78 | discussed in Section 2. These two sections are finalized by the description |
---|
79 | of the ``naive interpolation method'' that will be plugged into the coupler |
---|
80 | black box. This method have been coded in FORTRAN and is now under a |
---|
81 | testing stage. The results of the tests will be exposed in Part B of this |
---|
82 | letter. In Section 3, some steps towards the design of a less naive |
---|
83 | grid-to-grid interpolation method are indicated, and could be easily |
---|
84 | implemented in the ``naive method'' FORTRAN software, if they are |
---|
85 | retained as pertinent developments worth the effort. |
---|
86 | |
---|
87 | |
---|
88 | A lot of materials exposed here, at a pedagogical and |
---|
89 | simple level, are presented in a more complete and detailled manner in [3]. |
---|
90 | This is the case, for instance, for the gaussian basis and overlapping of |
---|
91 | Sections 2 and 3, and of the Lagrangian formalism. For furthers details, and |
---|
92 | in particular numerical tests of various interpolation methods in the 1D |
---|
93 | case, the reader is referred to [3]. |
---|
94 | |
---|
95 | |
---|
96 | |
---|
97 | |
---|
98 | |
---|
99 | |
---|
100 | \beginsection 2. BASIC FACTS ABOUT INTERPOLATION |
---|
101 | |
---|
102 | Interpolation is a problem that often occurs in physics, and we all have |
---|
103 | had some experience in it, in a way or another. Due to the numerous |
---|
104 | approachs that can be found in the literature, I have chosen to build a |
---|
105 | pedagogical presentation from scratch, restricted to the materials |
---|
106 | necessary to expose the ``naive method'' that I have coded. Simple |
---|
107 | notations for the indexing of masked grids are given, in order to be used |
---|
108 | in this presentation. The family of the linear interpolations associated to a |
---|
109 | grid is presented in a formalism focusing their functions basis set. |
---|
110 | Several illustrative examples are presented, ending with exposition of the |
---|
111 | ``$L$-neighbor naive method''. |
---|
112 | |
---|
113 | |
---|
114 | |
---|
115 | \beginsection 2.1 Grids and masks |
---|
116 | |
---|
117 | Let $r_{lm}= (x _{lm}, y _{lm})$, $ \; l=1,N_1$, $ \; m=1,N_2$ the |
---|
118 | coordinates of a rectangular 2D grid. In order to define the sea and the |
---|
119 | land, let $M _{lm}$ a $N_1 \times N_2$ matrix called ``mask'' such that, |
---|
120 | for instance, $M _{lm}=1$ means that $r_{lm}$ is a ``sea'' grid and $M |
---|
121 | _{lm}=0$ means that $r_{lm}$ is a ``land'' grid. In all the present |
---|
122 | discussion, we will only consider the sea points of the grid and discard the |
---|
123 | land points. If $N$ is the number of sea grid points we will denote them |
---|
124 | by $r_i=(x_i, y_i),\; i=1,N$ and no mask needs to be mentioned anymore in |
---|
125 | the presentation of the sea grid. The notion of mask is still useful |
---|
126 | however in the numerical implementation of algorithms dealing with |
---|
127 | thoses grids. Another advantage is the convenient generalization of the |
---|
128 | notations to 1D or 3D grids. |
---|
129 | |
---|
130 | |
---|
131 | |
---|
132 | |
---|
133 | |
---|
134 | \beginsection 2.2 The linear interpolation family of a grid |
---|
135 | |
---|
136 | The most common meaning of the word interpolation when dealing with a |
---|
137 | grid $(r_i),\; i=1,N$ with a real value $(a_i),\; i=1,N$ |
---|
138 | at each point is the |
---|
139 | determination of a function $f(r)$ on the whole space (1D, 2D or 3D |
---|
140 | depending on the case), which allows to extrapolate the values $a_i$'s in |
---|
141 | any point $r$ located on or outside the grid. |
---|
142 | Note that an interpolation does not necessarily requires $f(r_i)=a_i$. |
---|
143 | |
---|
144 | |
---|
145 | |
---|
146 | Any linear interpolation method can be expressed as |
---|
147 | $$ |
---|
148 | f(r)= \sum_{i=1}^{N} a_i \varphi_i(r) |
---|
149 | \eqno(2.1) |
---|
150 | $$ |
---|
151 | where the $ \varphi_i(r)$'s are functions localized around $r_i$ such that, |
---|
152 | for all $r$: |
---|
153 | $$ |
---|
154 | \sum_{i=1}^{N}\varphi_i(r) dr = 1 |
---|
155 | \eqno(2.2) |
---|
156 | $$ |
---|
157 | for all $i$. |
---|
158 | Let us call call these functions ``the basis'' of the interpolation. |
---|
159 | |
---|
160 | |
---|
161 | |
---|
162 | For a given grid, the family of all the possible interpolations of the form |
---|
163 | (2.1) is as big as all the possible choices of the basis functions. |
---|
164 | The choice of a particular interpolation depends on physical |
---|
165 | considerations, which depend on the use that one intends to make |
---|
166 | of the function $f(r)$. |
---|
167 | |
---|
168 | |
---|
169 | |
---|
170 | \medskip |
---|
171 | |
---|
172 | Let us assume that the |
---|
173 | general shape of the $ \varphi_i(r)$'s is more or less the same on all |
---|
174 | the grid. This is the case, for instance, when |
---|
175 | a given $ \varphi_i(r)$ ``covers'' |
---|
176 | roughly the same number of |
---|
177 | neighbors for most of the points of the grid. |
---|
178 | It is then possible to define the ``overlap order'' of an interpolation |
---|
179 | by looking at |
---|
180 | how much the support of the $ \varphi_i(r)$'s overlap, taking into account |
---|
181 | how fast the absolute value of $ \varphi_i(r)$. |
---|
182 | |
---|
183 | There is no need here to |
---|
184 | give a more precise definition that this intuitive concept of ``overlap |
---|
185 | order''. |
---|
186 | |
---|
187 | \vfill\eject |
---|
188 | |
---|
189 | |
---|
190 | \beginsection 2.3 The ``closest neighbor'' interpolation |
---|
191 | |
---|
192 | Given a distance in the space where stands the grid, for instance |
---|
193 | on a plane or on a |
---|
194 | sphere), the ``closest neighbor'' interpolation is such that |
---|
195 | $f(r)=a_{i(r)}$, where |
---|
196 | $a_{i(r)}$ is the value data at the closest neighbor |
---|
197 | $r_{i(r)}$, i.e. the grid point which minimizes the |
---|
198 | distances $(|r-r_i|), \; i=1,N$. |
---|
199 | |
---|
200 | |
---|
201 | |
---|
202 | This interpolation can be expressed under the form (2.1) with |
---|
203 | $\varphi_i(r)$ being 1 on the open domain $S_i$ of all the points |
---|
204 | $r$ which admit $r_i$ as their closest neighbor and $|S_i|$ its area. It |
---|
205 | is 1/2 on the boundary of this domain for the points at equal distance |
---|
206 | from 2 grid points, |
---|
207 | 1/4 at the ``corners'' at equal distance of 4 grid points, and so on. |
---|
208 | It is zero outside of the domain $S_i$. |
---|
209 | |
---|
210 | |
---|
211 | |
---|
212 | The overlap order of this interpolation is zero, because the support domains |
---|
213 | $S_i$ of the basis functions do not intersect. The function $f(r)$ of this |
---|
214 | interpolation is discontinuous. |
---|
215 | |
---|
216 | \vfill\eject |
---|
217 | |
---|
218 | \beginsection 2.4 Familiar Interpolations |
---|
219 | |
---|
220 | The linear interpolation is simple to understand on a 1D grid. Given $N$ |
---|
221 | values $(a_i) i=1,N$ at the real points $(x_i),\; i=1,N$ the linear |
---|
222 | interpolation reads: |
---|
223 | $$ |
---|
224 | f(x)=(a_{i+1}-a_i) { x-x_i \over x_{i+1} - x_i } + a_i |
---|
225 | \qquad \hbox{if} \quad x\in [x_i, x_{i+1}] \;, |
---|
226 | \eqno(2.3) |
---|
227 | $$ |
---|
228 | i.e. a piece linear function such that $f(x_i)=a_i$ for all points. |
---|
229 | |
---|
230 | |
---|
231 | |
---|
232 | Simple algebraic manipulations show that this interpolation can be put |
---|
233 | under the form (2.1) with |
---|
234 | $$ |
---|
235 | \varphi_i(x)=\left\{\matrix{ |
---|
236 | {x-x_{i-1}\over x_{i}-x_{i-1} } |
---|
237 | \qquad \hbox{if} \quad x\in [x_{i-1}, x_{i}] \cr |
---|
238 | {x-x_{i+1}\over x_{i}-x_{i+1} } |
---|
239 | \qquad \hbox{if} \quad x\in [x_{i}, x_{i+1}] \cr |
---|
240 | 0 |
---|
241 | \qquad \hbox{elsewhere}} \right. |
---|
242 | \eqno(2.4) |
---|
243 | $$ |
---|
244 | This function is equal to one at $x_i$, zero at the point $x_{i-1}$ and |
---|
245 | $x_{i+1}$ and outside $ [x_{i-1}, x_{i+1}]$, and is piecewise linear on two |
---|
246 | intervals. |
---|
247 | |
---|
248 | |
---|
249 | |
---|
250 | |
---|
251 | \bigskip |
---|
252 | |
---|
253 | The linear interpolation on a 2D grid $r_{lm}=(x_{lm},y_{lm})$ can also be |
---|
254 | put under the form |
---|
255 | $$ |
---|
256 | f(r)= \sum_{l,m} a_{lm} \varphi_{lm}(r) |
---|
257 | \eqno(2.5)\;, |
---|
258 | $$ |
---|
259 | with $\varphi_{lm}(r)$ equal to one at $r_{lm}$, zero at the four points |
---|
260 | $r_{l-1,m-1}$, $r_{l+1,m-1}$, $r_{l+1,m+1}$, $r_{l-1,m+1}$ and outside |
---|
261 | the polygon that they form, and piecewise linear on four triangles that |
---|
262 | they make with $r_{lm}$. |
---|
263 | |
---|
264 | |
---|
265 | |
---|
266 | \bigskip |
---|
267 | |
---|
268 | It can be said that the overlap order of the linear interpolation is one, |
---|
269 | as two basis functions can overlap on one grid mesh. It is intuitive to say |
---|
270 | that the higher the overlap, the smoother the interpolation. |
---|
271 | I will not try to consider here interpolations like quadratic interpolations |
---|
272 | or splines, which certainly enter into the form (2.1). |
---|
273 | |
---|
274 | |
---|
275 | |
---|
276 | \beginsection 2.5 ``Unique weight function'' interpolation |
---|
277 | |
---|
278 | Instead I will study a very natural interpolation base on the choice of a |
---|
279 | ``weight function'' $\phi(u)$, which can be, as an example, a gaussian: |
---|
280 | $$ |
---|
281 | \phi(u)= e^{-{u^2 \over 2\sigma^2}} \;. |
---|
282 | \eqno(2.6) |
---|
283 | $$ |
---|
284 | The normalization constant in front of the exponential can be anything as |
---|
285 | will be seen below. |
---|
286 | |
---|
287 | |
---|
288 | |
---|
289 | With the choice of a ``weight function'' one can define the interpolation |
---|
290 | $$ |
---|
291 | f(r)= {1\over Z(r)} \sum_{i=1}^{N} a_i\; \phi(r-r_i)\;, |
---|
292 | \eqno(2.7) |
---|
293 | $$ |
---|
294 | where |
---|
295 | $$ |
---|
296 | Z (r)= \sum_{i=1}^{N} \phi(r-r_i)\;, |
---|
297 | \eqno(2.8) |
---|
298 | $$ |
---|
299 | to ensure (2.2). |
---|
300 | This function $Z(r)$ depends on the geometry of the grid and on the weight |
---|
301 | function $\phi$, and its graphic visualization might give some valuable |
---|
302 | information about the chosen method. The same remark holds for the |
---|
303 | individual $ \varphi_i(r) = \phi(r-r_i) /Z(r)$. |
---|
304 | |
---|
305 | |
---|
306 | \vfill\eject |
---|
307 | |
---|
308 | |
---|
309 | \beginsection 2.6 The ``$L$ neighbors naive interpolation'' |
---|
310 | |
---|
311 | The numerical implementation of the above ``unique weight'' interpolation |
---|
312 | is simple to code, but it costs of order $N$ operation for each point $r$ |
---|
313 | where the value of $f(r)$ is calculated. In order to save computing time |
---|
314 | this interpolation can be approximated or mimiced by what I call the |
---|
315 | ``naive method of order $L$'': |
---|
316 | $$ |
---|
317 | f(r)={1\over Z(r)} \sum_{m\in I(r)} a_i \; \phi(r-r_m) |
---|
318 | \eqno(2.9) |
---|
319 | $$ |
---|
320 | where $I(r)=\{ i_1(r), i_2(r), ..., i_L(r) \} $ is the set containing the |
---|
321 | indices of the $L$ closest neighbors of $r$. The function $Z(r)$ ensures, |
---|
322 | as usual, that the sum of the weights is one, and is here: |
---|
323 | $$ |
---|
324 | Z(r)= \sum_{m\in I(r)} \phi(r-r_m) |
---|
325 | \eqno(2.10) |
---|
326 | $$ |
---|
327 | This function should also be worth graphic plots for various choices of the |
---|
328 | integer $L$ and the variance of $\phi$ ($\sigma$ for the gaussian). |
---|
329 | |
---|
330 | |
---|
331 | |
---|
332 | \medskip |
---|
333 | |
---|
334 | This ``naive interpolation'' can also be expressed under the form (2.1), |
---|
335 | with |
---|
336 | $$ |
---|
337 | \varphi_i(r) =\left\{\matrix{ \phi(r-r_m)/Z(r) \qquad \hbox{if} \quad |
---|
338 | r \in S^L_j \cr |
---|
339 | 0 \qquad \qquad \hbox{elsewhere}} \right. |
---|
340 | \eqno(2.11) |
---|
341 | $$ |
---|
342 | where $S_j^L$ is the domain of all the points $r$ which admit the grid |
---|
343 | point $r_j$ in the list of the $L$ closest neighbors. The values at the |
---|
344 | boundary has to be adjusted as for the ``closest neighbor'' interpolation, |
---|
345 | which is the particular case $L=1$ and $\phi=1$. |
---|
346 | |
---|
347 | |
---|
348 | |
---|
349 | |
---|
350 | |
---|
351 | |
---|
352 | |
---|
353 | \beginsection 3. GRID-TO-GRID INTERPOLATION |
---|
354 | |
---|
355 | There is a big difference between the interpolation from a set of data on a |
---|
356 | grid to a function $f(r)$ defined everywhere, and a grid-to-grid |
---|
357 | interpolation. In a grid-to-grid interpolation, a set of data is given on a |
---|
358 | first grid, the ``source grid'', and the interpolation must generate from |
---|
359 | these a set of data on a new grid, the ``target grid''. This problems occurs |
---|
360 | in the ocean-atmosphere coupling problem where fluxes must be transmitted |
---|
361 | from the atmospheric grid to the oceanic grid, and sea surface temperature |
---|
362 | (SST) |
---|
363 | from the oceanic grid to the atmospheric grid. |
---|
364 | Of course, given an interpolation on the source grid, one can apply the |
---|
365 | interpolated function $f(r)$ at each point of the target grid to define a |
---|
366 | grid-to-grid interpolation. But this choice can be extended through a |
---|
367 | Lagrangian formalism where the grid-to-grid interpolation depends on an |
---|
368 | interpolation on the source grid and also an interpolation of the target |
---|
369 | grid. Indeed the grid-to-grid interpolation depends on what the physicist |
---|
370 | intends to do with the interpolated values, and this can generally be |
---|
371 | expressed through the definition of an interpolation on the target grid. |
---|
372 | These considerations will guide the choices of the parameters of our ``naive'' |
---|
373 | grid-to-grid interpolation (number of neighbors and variance). |
---|
374 | |
---|
375 | \vfill\eject |
---|
376 | |
---|
377 | |
---|
378 | \beginsection 3.1 Sea domains and Ocean-Atmosphere coupling |
---|
379 | |
---|
380 | The grid of an Atmospheric General Circulation Model (AGCM) is made of |
---|
381 | ``land points'' and ``sea points''. The grid of an Oceanic General Circulation |
---|
382 | Model (OGCM) is, by definition, only made of ``oceanic points''. The AGCM |
---|
383 | ``sea points'' and the OGCM ``oceanic points'' generally do not coincide, |
---|
384 | except when the same grids has been chosen. |
---|
385 | |
---|
386 | |
---|
387 | |
---|
388 | The fluxes of the AGCM are known at all the points, land or sea, and one of |
---|
389 | the |
---|
390 | task of the coupler is to interpolate them for the OGCM grid. |
---|
391 | This can be done |
---|
392 | either by using all the grid, with no distinction between land and sea |
---|
393 | points, or only taking into account the sea points. This last point is still |
---|
394 | under debate. There is no such a choice |
---|
395 | however when passing the sea surface temperature (SST) of the OGCM grid |
---|
396 | to the AGCM ``sea points''. |
---|
397 | |
---|
398 | |
---|
399 | |
---|
400 | In all these cases, it is possible to denote by by $(r^a_i),\; i=1,N^a$ the |
---|
401 | ``source grid'' on which are known the data and $(r^b_i),\; i=1,N^b$ the |
---|
402 | ``target grid'' on which the interpolated values must be calculated. The |
---|
403 | two domains of the two grids are likely to have a large intersection and |
---|
404 | sometimes a not so small differences (the complementary of the |
---|
405 | intersection in the reunion of the the two domains). The influence of the |
---|
406 | difference set of the two grids in a grid-to-grid interpolation will be |
---|
407 | investigated in Part B of this letter, through numerical tests. |
---|
408 | |
---|
409 | |
---|
410 | |
---|
411 | \beginsection 3.2 Lagrangian formalism |
---|
412 | |
---|
413 | A grid-to-grid interpolation always assumes, even if not said, the |
---|
414 | choice of an interpolation on the source grid as well as on the target |
---|
415 | grid. |
---|
416 | |
---|
417 | |
---|
418 | |
---|
419 | Let |
---|
420 | $$ |
---|
421 | f(r)=\sum_{i=1}^{N^a} a_i \varphi_i(r) |
---|
422 | \eqno(3.1) |
---|
423 | $$ |
---|
424 | a chosen interpolation on the source grid A defined by the points $(r^a_i),\; |
---|
425 | i=1,N^a$, at which the values $a_i$ are known. |
---|
426 | |
---|
427 | |
---|
428 | |
---|
429 | Let |
---|
430 | $$ |
---|
431 | g(r)=\sum_{j=1}^{N^b} b_i \psi_j(r) |
---|
432 | \eqno(3.2) |
---|
433 | $$ |
---|
434 | a chosen interpolation on the target grid B defined by the points $(r^b_i),\; |
---|
435 | j=1,N^b$, at which the values $b_i$ are unknown an to be calculated |
---|
436 | through a grid-to-grid interpolation. |
---|
437 | |
---|
438 | |
---|
439 | |
---|
440 | Such a grid-to-grid interpolation can be defined by the minimization of |
---|
441 | the Lagrangian: |
---|
442 | $$ |
---|
443 | L(b_1, b_2, ..., b_{N^b}) = \int [f(r)-g(r)]^2 dr \;. |
---|
444 | \eqno(3.3) |
---|
445 | $$ |
---|
446 | The integral is on all the domain, owing the fact that the $ |
---|
447 | \varphi_i(r)$'s can be zero outside their support. |
---|
448 | It is natural to minimize this Lagrangian in order to have $f(r)$ as close |
---|
449 | as possible of $g(r)$. |
---|
450 | The set of coefficients $(b_j),\; j=1,N^b$ which realizes this minimization |
---|
451 | is found by saying that all the partial derivatives $\partial L / \partial |
---|
452 | b_j$ vanish. Some straightforward algebraic manipulations show that this |
---|
453 | lead the $N^b$ linear equations: |
---|
454 | $$ |
---|
455 | 2 \sum_{j=1}^{N^b} < \psi_m | \psi_j > b_j |
---|
456 | =2 \sum_{i=1}^{N^a} < \psi_m | \varphi_i> a_i |
---|
457 | \qquad \hbox{for} \quad m=1,N^b |
---|
458 | \eqno(3.4) |
---|
459 | $$ |
---|
460 | where |
---|
461 | $$ |
---|
462 | < \psi_m |\psi_j > = \int \psi_m (r)\psi_j(r) dr =: V_{mj} |
---|
463 | \eqno(3.5) |
---|
464 | $$ |
---|
465 | and |
---|
466 | $$ |
---|
467 | < \psi_m | \varphi_i > = \int \psi_m (r)\varphi_i(r) P(r) dr =: U_{mi} |
---|
468 | \;. |
---|
469 | \eqno(3.6) |
---|
470 | $$ |
---|
471 | |
---|
472 | |
---|
473 | |
---|
474 | In order to find the $N^b$ values $b_j$'s, one has to invert the linear |
---|
475 | system $VB = UA$: |
---|
476 | $$ |
---|
477 | \sum_{j=1}^{N_b} V_{ml}\; b_j = \sum_{i=1}^{N_a} U_{mi}\; a_i |
---|
478 | \qquad\hbox{for}\quad m=1,N^b |
---|
479 | \eqno(3.7) |
---|
480 | $$ |
---|
481 | |
---|
482 | |
---|
483 | |
---|
484 | The inversion of the $N_b$ x $N_b$ matrix $V$ is possible if the |
---|
485 | $\psi_j$'s are a good basis, and the grid-to-grid interpolation finally |
---|
486 | reads: $B=V^{-1}U A$. |
---|
487 | |
---|
488 | |
---|
489 | |
---|
490 | |
---|
491 | |
---|
492 | \beginsection 3.3 The choice of the two grid interpolations |
---|
493 | |
---|
494 | |
---|
495 | |
---|
496 | Choosing the ``closest neighbor'' interpolation for the target grid, in the |
---|
497 | above formalism, lead to a very easy implementation of the method because |
---|
498 | the matrix $V$ is equal to unity. The interpolation chosen on the source |
---|
499 | grid (grid A) is then the only one which determines the grid-to-grid |
---|
500 | interpolation and reads: |
---|
501 | $$ |
---|
502 | b_j= \sum _{i=1}^{N^a} a_i \overline{\varphi_i}^{S_j}(r_j^b) |
---|
503 | \eqno(3.8) |
---|
504 | $$ |
---|
505 | where |
---|
506 | $$ |
---|
507 | \overline{\varphi_i}^{S_j}(r_j^b)= \int_{S_j} \varphi_i(r) |
---|
508 | \eqno(3.9) |
---|
509 | $$ |
---|
510 | is the average of $\varphi$ on the ``closest neighbor'' cell surrounding |
---|
511 | $r_j$. When the variance of $\varphi_i$ is large compared to the size of |
---|
512 | this domain, this average is close to its value at $r_j$ and we have |
---|
513 | approximatively $b_j\sim f(r_j)$. On the contrary, if this variance is |
---|
514 | small compared to the size of the ``closest neighbor'' cells of the target |
---|
515 | grid (grid B), the grid-to-grid interpolation is approximatively the trivial |
---|
516 | interpolation for which all the $b_j$'s are equal to the mean value |
---|
517 | $\sum_i a_i$. |
---|
518 | |
---|
519 | |
---|
520 | |
---|
521 | |
---|
522 | \bigskip |
---|
523 | |
---|
524 | |
---|
525 | Having this in mind, I want to express here some intuitive ideas about |
---|
526 | the choice of the source and target grid interpolation, when building a |
---|
527 | grid-to-grid interpolation with the above Lagrangian formalism. Let us |
---|
528 | consider three situations. |
---|
529 | |
---|
530 | \medskip |
---|
531 | |
---|
532 | |
---|
533 | |
---|
534 | \item{*} {\bf Coarse to thin grid}. If the source grid A is coarser that the |
---|
535 | target grid B, it seems natural to choose the ``closest neighbor'' |
---|
536 | interpolation on grid B. On grid A, one feels that an smooth interpolation |
---|
537 | with an ``overlap order'' which is be 2 or 5 times |
---|
538 | to the thinest ratio of the |
---|
539 | two grids, i.e., the ratio of the average mesh size of the coarse grid A over |
---|
540 | the average mesh size of the thin grid B. |
---|
541 | |
---|
542 | \medskip |
---|
543 | |
---|
544 | \item{*} {\bf Thin to coarse grid}. It the source grid A is thiner that the |
---|
545 | target grid B. One can still try to choose the ``closest neighbor'' |
---|
546 | interpolation on grid B. |
---|
547 | On grid A, one feels that an smooth interpolation with an ``overlap order'' |
---|
548 | which is roughly equal to the coarsest ration of the two grid, i.e., the ratio |
---|
549 | of the average mesh size of the thin grid A over the average mesh size of |
---|
550 | the coarse grid B. |
---|
551 | |
---|
552 | |
---|
553 | |
---|
554 | |
---|
555 | |
---|
556 | \beginsection 3.4 The grid-to-grid ``naive'' method |
---|
557 | |
---|
558 | All the above considerations have been addressed in order to explain the |
---|
559 | grid-to-grid ``naive'' method that is being build as a temporary modulus, |
---|
560 | waiting for a less naive method to be coded. |
---|
561 | |
---|
562 | |
---|
563 | |
---|
564 | In this naive method the fluxes $b_j$ on the ocean grid B $(r^b_j),\; j=1,N^b$ |
---|
565 | are obtained from the fluxes $a_i$ of the atmospheric grid A $(r^a_i),\; |
---|
566 | i=1,N^a$ simply by applying the $L^a$ neighbors $\sigma_a$-gaussian |
---|
567 | interpolation with : |
---|
568 | $$ |
---|
569 | b_j={1\over Z^a(r^b_j)} \sum_{m\in I(r^b_j)} a_i e^{-{(r^b_j-r^a_m)^2 |
---|
570 | \over 2 \sigma_a^2} } |
---|
571 | \eqno(3.10) |
---|
572 | $$ |
---|
573 | where $I(r^b_j)=\{ i_1(r^b_j), i_2(r^b_j), ..., i_{L^a}(r^b_j) \} $ is the set |
---|
574 | containing the indices of $L^a$ closest neighbors A-grid points $r^a_m$ |
---|
575 | of the B-grid point $r^b_j$, and |
---|
576 | $$ |
---|
577 | Z^a(r^b_j) = \sum_{m\in I(r^b_j)} e^{-{(r^b_j-r^a_m)^2 \over 2 |
---|
578 | \sigma_a^2} }\;. |
---|
579 | \eqno(3.11) |
---|
580 | $$ |
---|
581 | |
---|
582 | |
---|
583 | |
---|
584 | \medskip |
---|
585 | |
---|
586 | When interpolating the SST (sea surface temperature) for the oceanic grid |
---|
587 | B to the atmospheric grid, the same principle is applied with |
---|
588 | the $L^b$ neighbors $\sigma_b$-gaussian interpolation: |
---|
589 | $$ |
---|
590 | a_i={1\over Z^b(r^a_i)} \sum_{m\in J(r^a_i)} b_i e^{-{(r^a_i-r^b_m)^2 |
---|
591 | \over 2 \sigma_b^2} } |
---|
592 | \eqno(3.12) |
---|
593 | $$ |
---|
594 | where $J(r^a_i)=\{ j_1(r^a_i), j_2(r^a_i), ..., i_{L^b}(r^a_i) \} $ is the set |
---|
595 | containing the indices of $L^b$ closest neighbors B-grid points $r^b_m$ |
---|
596 | of the A-grid point $r^a_i$, and |
---|
597 | $$ |
---|
598 | Z^b(r^a_i) = \sum_{m\in J(r^a_i)} e^{-{(r^a_i-r^b_m)^2 \over 2 |
---|
599 | \sigma_b^2} } |
---|
600 | \eqno(3.13) |
---|
601 | $$ |
---|
602 | |
---|
603 | \vfill\eject |
---|
604 | |
---|
605 | |
---|
606 | \bigskip |
---|
607 | |
---|
608 | The adjustable parameters of this ``naive method'' are the two integers |
---|
609 | $L^a$ and $L^b$ and the two variances $\sigma_a$ and $\sigma_b$. For a |
---|
610 | given couple of grids A and B, their choice can be guided by the above |
---|
611 | analysis base on the Lagrangian formalism. |
---|
612 | |
---|
613 | |
---|
614 | |
---|
615 | If Grid B is thiner than Grid A, which is the case for the |
---|
616 | Ocean-Atmosphere coupling, it is reasonable to choose $L^a$ or order 10 |
---|
617 | and $\sigma_a$ of the order of the average mesh of grid A. In the above |
---|
618 | formalism of the Lagrangian interpolation, this choice corresponds to |
---|
619 | the $L^a$ neighbors $\sigma_a$-gaussian interpolation on Grid A, i.e., |
---|
620 | and the ``delta-comb'' interpolation on Grid B, i.e., the choice of delta |
---|
621 | function for the $\psi_j$'s. The inverse Grid B to Grid A interpolation can |
---|
622 | be expressed the same way, but with a different choice of $L^b$. Indeed, it |
---|
623 | seems natural to choose $L^b$ such that $L^b/L^a$ is of the order of the |
---|
624 | ratio between the A-grid and B-grid average mesh size, and $\sigma_b$ of |
---|
625 | the order of $\sigma_a$. With this choice a composition between a A to B |
---|
626 | and B to A grid interpolation is close to the Identity, which is a reasonable |
---|
627 | criteria to impose when no other information is assumed on the nature of |
---|
628 | the fields to be interpolated. |
---|
629 | |
---|
630 | |
---|
631 | |
---|
632 | |
---|
633 | |
---|
634 | |
---|
635 | |
---|
636 | \beginsection 4. TOWARD A LESS NAIVE METHOD |
---|
637 | |
---|
638 | For the case of Ocean-Atmosphere coupling, information is known about |
---|
639 | the structure of the flux or SST fields to interpolated, and about the |
---|
640 | geometry of the grid. The naive method should be thus replaced by |
---|
641 | grid-to-grid interpolations adapted to the specific physical problem of the |
---|
642 | Ocean-Atmosphere coupling. These interpolation have to be tuned in order |
---|
643 | to answer specific quality criteria, imposed by physical arguments. This |
---|
644 | tuning can be achieved through the choice of the basis functions the |
---|
645 | source grid interpolation (the Atmospheric grid when passing fluxes, and |
---|
646 | the Oceanic one when passing the SST). Constraints like the conservation |
---|
647 | of the global integral of fluxes can be satisfied by adding Lagrange |
---|
648 | multipliers in the Lagrangian formalism. At last, a choice of the target |
---|
649 | grid interpolation, different from the ``Dirac comb'' interpolation can be |
---|
650 | necessary. In this last case the inversion of the matrix $V$ of Equation |
---|
651 | (3.6) has to be done once for all, and the multiplication by $V^{-1}$, which |
---|
652 | is generally not sparse, has to be performed for each new field to be |
---|
653 | interpolated. The numerical computation of the elements of $U$ by |
---|
654 | Monte-Carlo techniques is briefly discussed. |
---|
655 | |
---|
656 | |
---|
657 | |
---|
658 | \beginsection 4.2 Inhomogeneous and anisotropic interpolation |
---|
659 | |
---|
660 | The choice of the interpolation basis $\varphi^a_i(r)$ and $\varphi^b_j(r)$ |
---|
661 | of the grid A and B considered as source grids (i.e. when interpolating |
---|
662 | fluxes or SST respectively), can be tuned in function of the geometry of |
---|
663 | the two grids, and the knowledge of the physical structure of the fields to |
---|
664 | interpolate. For instance, some anisotropy could be given to these |
---|
665 | functions if the fields have a higher correlation in a particular direction |
---|
666 | (e.g. in latitudinal direction on the sphere). Near the edges of the two |
---|
667 | grids, i.e. the costs, a special choice of the basis function, adapted to the |
---|
668 | geometry or the physical processes should be performed. |
---|
669 | |
---|
670 | |
---|
671 | |
---|
672 | \beginsection 4.3 Lagrangian formalism with constraints |
---|
673 | |
---|
674 | When passing the flux of a quantity, like heat for instance, from the |
---|
675 | atmospheric grid to the oceanic grid it might be important to ensure that |
---|
676 | the amount of the heat lost by the AGCM is equal to the heat received by |
---|
677 | the OGCM. This conservation can be implemented by imposed one or |
---|
678 | several constraints in the definition of the interpolation. |
---|
679 | |
---|
680 | |
---|
681 | |
---|
682 | Taking the notation of the Langragian formalism of Section 3, the |
---|
683 | constraints |
---|
684 | $$ |
---|
685 | \int f(r) dr = \int g(r) dr |
---|
686 | \eqno(4.1) |
---|
687 | $$ |
---|
688 | defines a new grid-to-grid interpolation obtained by minimizing |
---|
689 | $$ |
---|
690 | L(b_1, b_2, ..., b_{N^b}) |
---|
691 | = \int [f(r)-g(r)]^2 dr + \lambda \int [f(r)-g(r)] dr |
---|
692 | \eqno(4.2) |
---|
693 | $$ |
---|
694 | where $\lambda$ is a Lagrangian multiplier. |
---|
695 | |
---|
696 | |
---|
697 | |
---|
698 | With this constraint, Equation (3.4) is now replaced by the system |
---|
699 | $$ |
---|
700 | \left\{ |
---|
701 | \eqalign{ |
---|
702 | 2 \sum_{j=1}^{N^b} < \psi_m | \psi_j > b_j |
---|
703 | =&2 \sum_{i=1}^{N^a} < \psi_m | \varphi_i> a_i - \lambda < \psi_m | 1 > |
---|
704 | \qquad \hbox{ for all } \quad m=1,N^b \cr |
---|
705 | \sum_{j=1}^{N^b} |
---|
706 | < \psi_j | 1 > b_j =& \sum_{i=1}^{N^a} < \varphi_i | 1 > |
---|
707 | a_i |
---|
708 | } \right. |
---|
709 | \eqno(4.3) |
---|
710 | $$ |
---|
711 | where |
---|
712 | $$ |
---|
713 | < \psi_j | 1 > = \int \psi_j dr = : H_j |
---|
714 | \eqno(4.4) |
---|
715 | $$ |
---|
716 | and |
---|
717 | $$ |
---|
718 | < \varphi_i | 1 > = \int \varphi_i dr = : G_j |
---|
719 | \eqno(4.4) |
---|
720 | $$ |
---|
721 | |
---|
722 | Equation (3.7) is now replaced by a system $VB=UA-\lambda H$ with $HB= |
---|
723 | GA$: |
---|
724 | $$ |
---|
725 | \left\{ |
---|
726 | \eqalign{ |
---|
727 | \sum_{j=1}^{N_b} V_{ml}\; b_j =& \sum_{i=1}^{N_a} U_{mi} \; a_i - \lambda |
---|
728 | G_m |
---|
729 | \qquad\hbox{for}\quad m=1,N^b \cr |
---|
730 | \sum_{j=1}^{N^b} H_j \; b_j =& \sum_{i=1}^{N^a} G_i\; a_i |
---|
731 | }\right. |
---|
732 | \eqno(4.5) |
---|
733 | $$ |
---|
734 | |
---|
735 | |
---|
736 | |
---|
737 | This system is solved by first inverting the matrix $V$, which leads to the |
---|
738 | system $B=V^{-1}UA-\lambda V^{-1}H$ with $HB= GA$, and by |
---|
739 | eliminating $\lambda$. Let us denote $B^*= V^{-1}UA$, which is the value |
---|
740 | of the coefficients that would be obtained without the constraint, and plug |
---|
741 | $B=B^*- \lambda V^{-1}H$ in the second equation. This gives the value of |
---|
742 | the Lagrange multiplier |
---|
743 | $$ |
---|
744 | \lambda= { H B^* - G A \over G V^{-1}H} \; |
---|
745 | \eqno(4.6) |
---|
746 | $$ |
---|
747 | which gives a measure of how strongly the constraint is violated by the first |
---|
748 | guess $B^*$. This violation error is redistributed in the solution of the |
---|
749 | problem: |
---|
750 | $$ |
---|
751 | \eqalign{ |
---|
752 | B=& B^* - {H B^* - G A \over H V^{-1}H} V^{-1}H \cr |
---|
753 | =& V^{-1}UA - {H V^{-1}UA - G A \over H V^{-1}H} V^{-1} H \;.} |
---|
754 | \eqno(4.7) |
---|
755 | $$ |
---|
756 | |
---|
757 | |
---|
758 | |
---|
759 | |
---|
760 | |
---|
761 | \beginsection 4.3 Practical considerations |
---|
762 | |
---|
763 | For the practical implementation of the interpolation methods which are |
---|
764 | derived from a Lagrangian formalism, one needs to know the numerical |
---|
765 | values of the quantities $V_{mj}=~~~<\psi_m | \psi_j> $, $U_{mi}= <\psi_m | |
---|
766 | \phi_i> $, $H_m= <\psi_m | 1> $ and $G_i= <\varphi_i | 1> $. Even for the |
---|
767 | next neighbor interpolation basis, the numerical evaluation of these |
---|
768 | numbers is not straightforward. The complexity of this problem, however, |
---|
769 | do not exceed the one of the determination of the double integral of an |
---|
770 | arbitrary function, which is an often encountered problem in physics. |
---|
771 | Various methods can be found depending of the particular situation. |
---|
772 | |
---|
773 | |
---|
774 | One of them is the evaluation of these integrals through Monte-Carlo |
---|
775 | techniques. For instance, if the functions $\psi_m$'s are the next neighbor |
---|
776 | interpolation basis, the $H_m$'s are the surfaces of the non-intersecting |
---|
777 | domains $S_m$'s (see Section 2.3). If the surfaces of the grid meshes are |
---|
778 | known, and if the grid is not too distorted (in a sense to be defined later), |
---|
779 | it is possible to reduce the number of random shootings. Inside each grid |
---|
780 | mesh formed by the four grid points $r_1$, $r_2$, $r_3$ and $r_4$, one |
---|
781 | can choose at random, with a uniform probability, a big number of points. |
---|
782 | For each of these points $r$, it is only necessary to calculate the four |
---|
783 | values $\psi_1(r)$, $\psi_2(r)$, $\psi_3(r)$ and $\psi_4(r)$, associated to |
---|
784 | the four grid points at the corners of the mesh, provided that the grid is |
---|
785 | not too much distorted (by definition). If so, the random point has to be |
---|
786 | chosen uniformly on a surface larger than the considered mesh. The same |
---|
787 | algorithm can be apply for the calculation the the $U_{mi}$'s. If the |
---|
788 | $\varphi_i$'s are the basis function of the $L$-closest neighbors |
---|
789 | interpolation, this algorithm can be generalized. |
---|
790 | |
---|
791 | \vfill\eject |
---|
792 | |
---|
793 | \beginsection 5. CONCLUSION |
---|
794 | |
---|
795 | The goal of this letter was to expose the ``naive'' grid-to-grid |
---|
796 | interpolation method that I have coded for the interpolation black box of |
---|
797 | the coupler. For each point of the target grid, the algorithm search for its |
---|
798 | $L$-closest neighbor and calculates weight coefficients through the |
---|
799 | evaluation of a |
---|
800 | gaussian function. This has to be done once for all, given a couple of grid, |
---|
801 | the addressed and the weight coefficients being the only information to be |
---|
802 | stored. If $L$ is not big, this storage is not huge, and the number of |
---|
803 | operation is also reasonable. Test of this methods with several type of |
---|
804 | grids will be presented in Part B. |
---|
805 | |
---|
806 | In fact, the rough exposition of the method could have been done quite |
---|
807 | more shortly if this had been the only considered goal. But the different |
---|
808 | formalisms presented in this letter are candidate to provide a framework |
---|
809 | to the intuitive ideas that one often generates when dealing with this |
---|
810 | interpolation problem. It can also be useful to try to translate into these |
---|
811 | formalisms more complicated interpolation methods, and, if not possible, |
---|
812 | to see which of their features do not enter in this framework. |
---|
813 | |
---|
814 | |
---|
815 | |
---|
816 | \bigskip |
---|
817 | |
---|
818 | \noindent {\bf Acknowledgments:} I thank all the participants of the ``Adm |
---|
819 | \& Science Team Meeting'' of March 11th, for their remarks which have |
---|
820 | useful for the drafting of this letter. |
---|
821 | |
---|
822 | |
---|
823 | |
---|
824 | \beginsection REFERENCES |
---|
825 | |
---|
826 | \def\ref{\parskip 12pt \leftskip 20pt \parindent -20pt} |
---|
827 | \def\endref{\parskip 0pt \leftskip 0pt \parindent 20pt} |
---|
828 | |
---|
829 | |
---|
830 | |
---|
831 | \ref |
---|
832 | [1] |
---|
833 | O. THUAL, Gathering information for a coupler, |
---|
834 | {\it Epicoa \ } {\bf 1009} (1991). |
---|
835 | |
---|
836 | |
---|
837 | |
---|
838 | [2] |
---|
839 | L. TERRAY, A simple communication scheme between coupler and GCM's, |
---|
840 | {\it Epicoa \ } {\bf 0203} (1992). |
---|
841 | |
---|
842 | |
---|
843 | |
---|
844 | [3] |
---|
845 | N. H. BARTH, A Conservative Scheme for Passing Variables Between |
---|
846 | Coupled Models of the Ocean an Atmosphere, |
---|
847 | {\it Technical Report\ }, CERFACS (1992). |
---|
848 | |
---|
849 | |
---|
850 | |
---|
851 | [4] ``Appel d'id\'ees pour la cr\'eation du coupleur du mod\`ele |
---|
852 | communautaire'', written by L. Terray, October 10th 1991. |
---|
853 | |
---|
854 | |
---|
855 | |
---|
856 | \endref |
---|
857 | |
---|
858 | |
---|
859 | |
---|
860 | \end |
---|
861 | |
---|
862 | |
---|
863 | |
---|
864 | |
---|
865 | |
---|