Context Navigation

compute_fromreg_imoms3_weigaddr.pro @ 163

Last change on this file since 163 was 163, checked in by navarro, 18 years ago
header improvements : type of parameters and keywords, default values, spell checking + idldoc assistant (IDL online_help)
Property svn:keywords set to `Id`
File size: 16.4 KB

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

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

Download in other formats:

Original Format