Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 138

Last change on this file since 138 was 136, checked in by pinsard, 18 years ago
some improvements and corrections in some .pro file according to aspell and idldoc log file
Property svn:keywords set to `Id`
File size: 16.4 KB

Line
1	;+
2	;
3	; @file_comments
4	; compute the weight and address neede to interpolate data from a
5	; "regular grid" to any grid using the imoms3 method
6	;
7	; @categories interpolation
8	;
9	; @param alonin {in}{required}
10	; longitude of the input data
11	;
12	; @param alatin {in}{required}
13	; latitude of the input data
14	;
15	; @param olonin {in}{required}
16	; longitude of the output data
17	; @param olat {in}{required}
18	; latitude of the output data
19	;
20	; @keyword NONORTHERNLINE
21	; @keyword NOSOUTHERNLINE
22	; activate if you don't want to take into account the northen/southern line
23	; of the input data when perfoming the interpolation.
24	;
25	; @param weig {out}
26	; @param addr {out}
27	; 2D arrays, weig and addr are the weight and addresses used to
28	; perform the interpolation:
29	; dataout = total(weig*datain[addr], 1)
30	; dataout = reform(dataout, jpio, jpjo, /over)
31	;
32	; @restrictions
33	; - the input grid must be a "regular/rectangular grid", defined as a grid for
34	; which each lontitudes lines have the same latitude and each latitudes columns
35	; have the same longitude.
36	; - We supposed the data are located on a sphere, with a periodicity along
37	; the longitude.
38	; - points located between the first/last 2 lines are interpolated
39	; using a imoms3 interpolation along the longitudinal direction and linear
40	; interpolation along the latitudinal direction
41	; - points located out of the southern and northern boundaries are interpolated
42	; using a imoms3 interpolation only along the longitudinal direction.
43	;
44	; @history
45	; November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)
46	; March 2006: works for rectangular grids
47	;
48	; @version $Id$
49	;
50	;-
51	;
52	;----------------------------------------------------------
53	;----------------------------------------------------------
54	;
55	PRO compute_fromreg_imoms3_weigaddr, alonin, alatin, olonin, olat, weig, addr $
56	, NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline
57	;
58	compile_opt idl2, strictarrsubs
59	;
60	alon = alonin
61	alat = alatin
62	olon = olonin
63	;
64	jpia = n_elements(alon)
65	jpja = n_elements(alat)
66	;
67	jpio = (size(olon, /dimensions))[0]
68	jpjo = (size(olon, /dimensions))[1]
69	;
70	; alon
71	minalon = min(alon, max = maxalon)
72	IF maxalon-minalon GE 360. THEN stop
73	; alon must be monotonically increasing
74	IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN
75	shiftx = -(where(alon EQ min(alon)))[0]
76	alon = shift(alon, shiftx)
77	IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop
78	ENDIF ELSE shiftx = 0
79	; alon is it regularly spaced?
80	step = alon-shift(alon, 1)
81	step[0] = step[0] + 360.
82	IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregx = 1
83	; we extend the longitude range of alon (-> easy interpolation even
84	; near minalon et maxalon)
85	toadd = 10*jpia/360+1
86	alon = [alon[jpia-toadd:jpia-1]-360., alon[*], alon[0:toadd-1]+360.]
87	jpia = jpia+2*toadd
88	; alat
89	revy = alat[0] GT alat[1]
90	IF revy THEN alat = reverse(alat)
91	; alat must be monotonically increasing
92	IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop
93	; alat is it regularly spaced?
94	step = alat-shift(alat, 1)
95	step = step[1:jpja - 1L]
96	IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregy = 1
97	;
98	if keyword_set(nonorthernline) then BEGIN
99	jpja = jpja - 1L
100	alat = alat[0: jpja-1L]
101	ENDIF
102	if keyword_set(nosouthernline) then BEGIN
103	alat = alat[1: jpja-1L]
104	jpja = jpja - 1L
105	ENDIF
106	; olon between minalon et minalon+360
107	out = where(olon LT minalon)
108	WHILE out[0] NE -1 DO BEGIN
109	olon[out] = olon[out]+360.
110	out = where(olon LT minalon)
111	ENDWHILE
112	out = where(olon GE minalon+360.)
113	WHILE out[0] NE -1 DO BEGIN
114	olon[out] = olon[out]- 360.
115	out = where(olon GE minalon+360.)
116	ENDWHILE
117	; make sure that all values of olon are located within values of alon
118	IF min(olon, max = ma) LT minalon THEN stop
119	IF ma GE minalon+360. THEN stop
120	;
121	xaddr = lonarr(16, jpio*jpjo)
122	yaddr = lonarr(16, jpio*jpjo)
123	weig = fltarr(16, jpio*jpjo)
124	;
125	indexlon = value_locate(alon, olon)
126	IF total(alon[indexlon] GT olon) NE 0 THEN stop
127	IF total(alon[indexlon + 1L] LE olon) NE 0 THEN stop
128	IF (where(indexlon LE 1L ))[0] NE -1 THEN stop
129	IF (where(indexlon GE jpia-3L))[0] NE -1 THEN stop
130	indexlat = value_locate(alat, olat)
131	;
132	; for the ocean points located below the atm line
133	; jpja-2 and above the line 1
134	; for those points we can always find 16 neighbors
135	; imoms interpolation along longitude and latitude
136	;
137	short = where(indexlat LT jpja-2L AND indexlat GE 1L)
138	ilon = indexlon[short]
139	ilat = indexlat[short]
140	;
141	IF NOT keyword_set(noregy) THEN BEGIN
142	delta = alat[ilat+1L]-alat[ilat]
143	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
144	delta = delta[0]
145	;
146	d0 = (alat[ilat-1L]-olat[short])/delta
147	IF min(d0, max = ma) LE -2 THEN stop
148	IF ma GT -1 THEN stop
149	wy0 = imoms3(temporary(d0))
150	d1 = (alat[ilat ]-olat[short])/delta
151	IF min(d1, max = ma) LE -1 THEN stop
152	IF ma GT 0 THEN stop
153	wy1 = imoms3(temporary(d1))
154	d2 = (alat[ilat+1L]-olat[short])/delta
155	IF min(d2, max = ma) LE 0 THEN stop
156	IF ma GT 1 THEN stop
157	wy2 = imoms3(temporary(d2))
158	d3 = (alat[ilat+2L]-olat[short])/delta
159	IF min(d3, max = ma) LE 1 THEN stop
160	IF ma GT 2 THEN stop
161	wy3 = imoms3(temporary(d3))
162	ENDIF ELSE BEGIN
163	nele = n_elements(short)
164	wy0 = fltarr(nele)
165	wy1 = fltarr(nele)
166	wy2 = fltarr(nele)
167	wy3 = fltarr(nele)
168	FOR i = 0L, nele-1 DO BEGIN
169	IF i MOD 10000 EQ 0 THEN print, i
170	newlat = spl_incr(alat[ilat[i]-1L:ilat[i]+2L], [-1., 0., 1., 2.], olat[short[i]])
171	IF newlat LE 0 THEN stop
172	IF newlat GT 1 THEN stop
173	wy0[i] = imoms3(newlat+1)
174	wy1[i] = imoms3(newlat)
175	wy2[i] = imoms3(1-newlat)
176	wy3[i] = imoms3(2-newlat)
177	ENDFOR
178	ENDELSE
179	;
180	mi = min(wy0+wy1+wy2+wy3, max = ma)
181	IF abs(mi-1) GE 1.e-6 THEN stop
182	IF abs(ma-1) GE 1.e-6 THEN stop
183	;
184	IF NOT keyword_set(noregx) THEN BEGIN
185	delta = alon[ilon]-alon[ilon-1L]
186	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
187	delta = delta[0]
188	;
189	d0 = (alon[ilon-1L]-olon[short])/delta
190	IF min(d0, max = ma) LE -2 THEN stop
191	IF ma GT -1 THEN stop
192	wx0 = imoms3(temporary(d0))
193	d1 = (alon[ilon ]-olon[short])/delta
194	IF min(d1, max = ma) LE -1 THEN stop
195	IF ma GT 0 THEN stop
196	wx1 = imoms3(temporary(d1))
197	d2 = (alon[ilon+1L]-olon[short])/delta
198	IF min(d2, max = ma) LE 0 THEN stop
199	IF ma GT 1 THEN stop
200	wx2 = imoms3(temporary(d2))
201	d3 = (alon[ilon+2L]-olon[short])/delta
202	IF min(d3, max = ma) LE 1 THEN stop
203	IF ma GT 2 THEN stop
204	wx3 = imoms3(temporary(d3))
205	ENDIF ELSE BEGIN
206	nele = n_elements(short)
207	wx0 = fltarr(nele)
208	wx1 = fltarr(nele)
209	wx2 = fltarr(nele)
210	wx3 = fltarr(nele)
211	FOR i = 0L, nele-1 DO BEGIN
212	IF i MOD 10000 EQ 0 THEN print, i
213	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
214	IF newlon LE 0 THEN stop
215	IF newlon GT 1 THEN stop
216	wx0[i] = imoms3(newlon+1)
217	wx1[i] = imoms3(newlon)
218	wx2[i] = imoms3(1-newlon)
219	wx3[i] = imoms3(2-newlon)
220	ENDFOR
221	ENDELSE
222	;
223	mi = min(wx0+wx1+wx2+wx3, max = ma)
224	IF abs(mi-1) GE 1.e-6 THEN stop
225	IF abs(ma-1) GE 1.e-6 THEN stop
226	;
227	; line 0
228	xaddr[0, short] = ilon - 1L
229	xaddr[1, short] = ilon
230	xaddr[2, short] = ilon + 1L
231	xaddr[3, short] = ilon + 2L
232	yaddr[0, short] = ilat - 1L
233	yaddr[1, short] = yaddr[0, short]
234	yaddr[2, short] = yaddr[0, short]
235	yaddr[3, short] = yaddr[0, short]
236	weig[0, short] = wx0 * wy0
237	weig[1, short] = wx1 * wy0
238	weig[2, short] = wx2 * wy0
239	weig[3, short] = wx3 * wy0
240	; line 1
241	xaddr[4, short] = ilon - 1L
242	xaddr[5, short] = ilon
243	xaddr[6, short] = ilon + 1L
244	xaddr[7, short] = ilon + 2L
245	yaddr[4, short] = ilat
246	yaddr[5, short] = ilat
247	yaddr[6, short] = ilat
248	yaddr[7, short] = ilat
249	weig[4, short] = wx0 * wy1
250	weig[5, short] = wx1 * wy1
251	weig[6, short] = wx2 * wy1
252	weig[7, short] = wx3 * wy1
253	; line 2
254	xaddr[8, short] = ilon - 1L
255	xaddr[9, short] = ilon
256	xaddr[10, short] = ilon + 1L
257	xaddr[11, short] = ilon + 2L
258	yaddr[8, short] = ilat + 1L
259	yaddr[9, short] = yaddr[8, short]
260	yaddr[10, short] = yaddr[8, short]
261	yaddr[11, short] = yaddr[8, short]
262	weig[8, short] = wx0 * wy2
263	weig[9, short] = wx1 * wy2
264	weig[10, short] = wx2 * wy2
265	weig[11, short] = wx3 * wy2
266	; line 3
267	xaddr[12, short] = ilon - 1L
268	xaddr[13, short] = ilon
269	xaddr[14, short] = ilon + 1L
270	xaddr[15, short] = ilon + 2L
271	yaddr[12, short] = ilat + 2L
272	yaddr[13, short] = yaddr[12, short]
273	yaddr[14, short] = yaddr[12, short]
274	yaddr[15, short] = yaddr[12, short]
275	weig[12, short] = wx0 * wy3
276	weig[13, short] = wx1 * wy3
277	weig[14, short] = wx2 * wy3
278	weig[15, short] = wx3 * wy3
279	;
280	mi = min(total(weig[*, short], 1), max = ma)
281	IF abs(mi-1) GE 1.e-6 THEN stop
282	IF abs(ma-1) GE 1.e-6 THEN stop
283	;
284	; for the ocean points located between the atm lines
285	; jpja-2 and jpja-1 or between the atm lines 0 and 1
286	; linear interpolation between line 1 and line 2
287	;
288	short = where(indexlat EQ jpja-2L OR indexlat EQ 0)
289	IF short[0] NE -1 THEN BEGIN
290	ilon = indexlon[short]
291	ilat = indexlat[short]
292	;
293	delta = alat[ilat+1L]-alat[ilat]
294	IF NOT keyword_set(noregy) THEN BEGIN
295	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
296	delta = delta[0]
297	ENDIF
298	;
299	d1 = (alat[ilat ]-olat[short])/delta
300	IF min(d1, max = ma) LE -1 THEN stop
301	IF ma GT 0 THEN stop
302	wy1 = 1.+ temporary(d1)
303	d2 = (alat[ilat+1L]-olat[short])/delta
304	IF min(d2, max = ma) LE 0 THEN stop
305	IF ma GT 1 THEN stop
306	wy2 = 1.- temporary(d2)
307	;
308	mi = min(wy1+wy2, max = ma)
309	IF abs(mi-1) GE 1.e-6 THEN stop
310	IF abs(ma-1) GE 1.e-6 THEN stop
311	; but imoms3 along the longitude
312	IF NOT keyword_set(noregx) THEN BEGIN
313	delta = alon[ilon]-alon[ilon-1L]
314	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
315	delta = delta[0]
316	;
317	d0 = (alon[ilon-1L]-olon[short])/delta
318	IF min(d0, max = ma) LE -2 THEN stop
319	IF ma GT -1 THEN stop
320	wx0 = imoms3(temporary(d0))
321	d1 = (alon[ilon ]-olon[short])/delta
322	IF min(d1, max = ma) LE -1 THEN stop
323	IF ma GT 0 THEN stop
324	wx1 = imoms3(temporary(d1))
325	d2 = (alon[ilon+1L]-olon[short])/delta
326	IF min(d2, max = ma) LE 0 THEN stop
327	IF ma GT 1 THEN stop
328	wx2 = imoms3(temporary(d2))
329	d3 = (alon[ilon+2L]-olon[short])/delta
330	IF min(d3, max = ma) LE 1 THEN stop
331	IF ma GT 2 THEN stop
332	wx3 = imoms3(temporary(d3))
333	ENDIF ELSE BEGIN
334	nele = n_elements(short)
335	wx0 = fltarr(nele)
336	wx1 = fltarr(nele)
337	wx2 = fltarr(nele)
338	wx3 = fltarr(nele)
339	FOR i = 0L, nele-1 DO BEGIN
340	IF i MOD 10000 EQ 0 THEN print, i
341	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
342	IF newlon LE 0 THEN stop
343	IF newlon GT 1 THEN stop
344	wx0[i] = imoms3(newlon+1)
345	wx1[i] = imoms3(newlon)
346	wx2[i] = imoms3(1-newlon)
347	wx3[i] = imoms3(2-newlon)
348	ENDFOR
349	ENDELSE
350	;
351	mi = min(wx0+wx1+wx2+wx3, max = ma)
352	IF abs(mi-1) GE 1.e-6 THEN stop
353	IF abs(ma-1) GE 1.e-6 THEN stop
354	; line 1
355	xaddr[0, short] = ilon - 1L
356	xaddr[1, short] = ilon
357	xaddr[2, short] = ilon + 1L
358	xaddr[3, short] = ilon + 2L
359	yaddr[0, short] = ilat
360	yaddr[1, short] = ilat
361	yaddr[2, short] = ilat
362	yaddr[3, short] = ilat
363	weig[0, short] = wx0 * wy1
364	weig[1, short] = wx1 * wy1
365	weig[2, short] = wx2 * wy1
366	weig[3, short] = wx3 * wy1
367	; line 2
368	xaddr[4, short] = ilon - 1L
369	xaddr[5, short] = ilon
370	xaddr[6, short] = ilon + 1L
371	xaddr[7, short] = ilon + 2L
372	yaddr[4, short] = ilat + 1L
373	yaddr[5, short] = yaddr[4, short]
374	yaddr[6, short] = yaddr[4, short]
375	yaddr[7, short] = yaddr[4, short]
376	weig[4, short] = wx0 * wy2
377	weig[5, short] = wx1 * wy2
378	weig[6, short] = wx2 * wy2
379	weig[7, short] = wx3 * wy2
380	;
381	mi = min(total(weig[*, short], 1), max = ma)
382	IF abs(mi-1) GE 1.e-6 THEN stop
383	IF abs(ma-1) GE 1.e-6 THEN stop
384	;
385	ENDIF
386	;
387	; for the ocean points located below the line 0
388	; Interpolation only along the longitude
389	;
390	short = where(indexlat EQ -1)
391	IF short[0] NE -1 THEN BEGIN
392	ilon = indexlon[short]
393	;
394	IF NOT keyword_set(noregx) THEN BEGIN
395	delta = alon[ilon]-alon[ilon-1L]
396	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
397	delta = delta[0]
398	;
399	d0 = (alon[ilon-1L]-olon[short])/delta
400	IF min(d0, max = ma) LE -2 THEN stop
401	IF ma GT -1 THEN stop
402	wx0 = imoms3(temporary(d0))
403	d1 = (alon[ilon ]-olon[short])/delta
404	IF min(d1, max = ma) LE -1 THEN stop
405	IF ma GT 0 THEN stop
406	wx1 = imoms3(temporary(d1))
407	d2 = (alon[ilon+1L]-olon[short])/delta
408	IF min(d2, max = ma) LE 0 THEN stop
409	IF ma GT 1 THEN stop
410	wx2 = imoms3(temporary(d2))
411	d3 = (alon[ilon+2L]-olon[short])/delta
412	IF min(d3, max = ma) LE 1 THEN stop
413	IF ma GT 2 THEN stop
414	wx3 = imoms3(temporary(d3))
415	ENDIF ELSE BEGIN
416	nele = n_elements(short)
417	wx0 = fltarr(nele)
418	wx1 = fltarr(nele)
419	wx2 = fltarr(nele)
420	wx3 = fltarr(nele)
421	FOR i = 0L, nele-1 DO BEGIN
422	IF i MOD 10000 EQ 0 THEN print, i
423	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
424	IF newlon LE 0 THEN stop
425	IF newlon GT 1 THEN stop
426	wx0[i] = imoms3(newlon+1)
427	wx1[i] = imoms3(newlon)
428	wx2[i] = imoms3(1-newlon)
429	wx3[i] = imoms3(2-newlon)
430	ENDFOR
431	ENDELSE
432	;
433	mi = min(wx0+wx1+wx2+wx3, max = ma)
434	IF abs(mi-1) GE 1.e-6 THEN stop
435	IF abs(ma-1) GE 1.e-6 THEN stop
436	; line 1
437	xaddr[0, short] = ilon - 1L
438	xaddr[1, short] = ilon
439	xaddr[2, short] = ilon + 1L
440	xaddr[3, short] = ilon + 2L
441	yaddr[0:3, short] = 0
442	weig[0, short] = wx0
443	weig[1, short] = wx1
444	weig[2, short] = wx2
445	weig[3, short] = wx3
446	;
447	mi = min(total(weig[*, short], 1), max = ma)
448	IF abs(mi-1) GE 1.e-6 THEN stop
449	IF abs(ma-1) GE 1.e-6 THEN stop
450	;
451	ENDIF
452	;
453	; for the ocean points located above jpia-1
454	; Interpolation only along the longitude
455	;
456	short = where(indexlat EQ jpja-1L)
457	IF short[0] NE -1 THEN BEGIN
458	ilon = indexlon[short]
459	;
460	IF NOT keyword_set(noregx) THEN BEGIN
461	delta = alon[ilon]-alon[ilon-1L]
462	IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop
463	delta = delta[0]
464	;
465	d0 = (alon[ilon-1L]-olon[short])/delta
466	IF min(d0, max = ma) LE -2 THEN stop
467	IF ma GT -1 THEN stop
468	wx0 = imoms3(temporary(d0))
469	d1 = (alon[ilon ]-olon[short])/delta
470	IF min(d1, max = ma) LE -1 THEN stop
471	IF ma GT 0 THEN stop
472	wx1 = imoms3(temporary(d1))
473	d2 = (alon[ilon+1L]-olon[short])/delta
474	IF min(d2, max = ma) LE 0 THEN stop
475	IF ma GT 1 THEN stop
476	wx2 = imoms3(temporary(d2))
477	d3 = (alon[ilon+2L]-olon[short])/delta
478	IF min(d3, max = ma) LE 1 THEN stop
479	IF ma GT 2 THEN stop
480	wx3 = imoms3(temporary(d3))
481	ENDIF ELSE BEGIN
482	nele = n_elements(short)
483	wx0 = fltarr(nele)
484	wx1 = fltarr(nele)
485	wx2 = fltarr(nele)
486	wx3 = fltarr(nele)
487	FOR i = 0L, nele-1 DO BEGIN
488	IF i MOD 10000 EQ 0 THEN print, i
489	newlon = spl_incr(alon[ilon[i]-1L:ilon[i]+2L], [-1., 0., 1., 2.], olon[short[i]])
490	IF newlon LE 0 THEN stop
491	IF newlon GT 1 THEN stop
492	wx0[i] = imoms3(newlon+1)
493	wx1[i] = imoms3(newlon)
494	wx2[i] = imoms3(1-newlon)
495	wx3[i] = imoms3(2-newlon)
496	ENDFOR
497	ENDELSE
498	;
499	mi = min(wx0+wx1+wx2+wx3, max = ma)
500	IF abs(mi-1) GE 1.e-6 THEN stop
501	IF abs(ma-1) GE 1.e-6 THEN stop
502	; line 1
503	xaddr[0, short] = ilon-1L
504	xaddr[1, short] = ilon
505	xaddr[2, short] = ilon+1L
506	xaddr[3, short] = ilon+2L
507	yaddr[0:3, short] = jpja-1L
508	weig[0, short] = wx0
509	weig[1, short] = wx1
510	weig[2, short] = wx2
511	weig[3, short] = wx3
512	;
513	mi = min(total(weig[*, short], 1), max = ma)
514	IF abs(mi-1) GE 1.e-6 THEN stop
515	IF abs(ma-1) GE 1.e-6 THEN stop
516	;
517	ENDIF
518	;
519	; Come back to the original index of atm grid without longitudinal overlap.
520	;
521	;
522	xaddr = temporary(xaddr) - toadd
523	jpia = jpia - 2*toadd
524	; make sure all values are ge 0
525	xaddr = temporary(xaddr) + jpia
526	; range the values between 0 and jpia-1
527	xaddr = temporary(xaddr) mod jpia
528	;
529	; take into account shiftx if needed
530	IF shiftx NE 0 THEN xaddr = (temporary(xaddr) - shiftx) MOD jpia
531	; take into account nosouthernline and nonorthernline
532	if keyword_set(nosouthernline) then BEGIN
533	yaddr = temporary(yaddr) + 1L
534	jpja = jpja + 1L
535	ENDIF
536	if keyword_set(nonorthernline) then jpja = jpja + 1L
537	; take into account revy if needed
538	IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr)
539	; ;
540	addr = temporary(yaddr)*jpia+temporary(xaddr)
541	;
542	RETURN
543	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format