1 | ;+ |
---|
2 | ; NAME: interp.pro |
---|
3 | ; |
---|
4 | ; PURPOSE: Performs a bilinear or bicubic interpolation on |
---|
5 | ; data from a regular origin grid to any grid |
---|
6 | |
---|
7 | ; CATEGORY: Subroutine |
---|
8 | ; |
---|
9 | ; CALLING SEQUENCE: interp, znumdta, zdata_name, mask, t, glam, gphi, zresul |
---|
10 | ; |
---|
11 | ; INPUTS: |
---|
12 | ; |
---|
13 | ; znumdta : ID of input data file |
---|
14 | ; zdata_name : input data name |
---|
15 | ; mask : data mask |
---|
16 | ; t : timestep |
---|
17 | ; glam : longitude of target grid |
---|
18 | ; gphi : latitude " " " |
---|
19 | ; |
---|
20 | ; KEYWORD PARAMETERS: None |
---|
21 | ; |
---|
22 | ; OUTPUTS: |
---|
23 | ; zresul : interpolated field |
---|
24 | ; |
---|
25 | ; COMMON BLOCKS: |
---|
26 | ; common_interp.pro |
---|
27 | ; |
---|
28 | ; SIDE EFFECTS: |
---|
29 | ; |
---|
30 | ; RESTRICTIONS: |
---|
31 | ; |
---|
32 | ; EXAMPLE: |
---|
33 | ; |
---|
34 | ; MODIFICATION HISTORY: 11/99 A. Jouzeau |
---|
35 | ; 04/00 A. Jouzeau, M. Imbard |
---|
36 | ; 08/00 A. Jouzeau |
---|
37 | ; 10/01 R. Hordoir : Added boundary conditions, |
---|
38 | ; Bug correction, North Pole Problem |
---|
39 | ; 02/02 R. Hordoir : Added North Pole Treatment |
---|
40 | ; for non coherent winds. |
---|
41 | ; 03/2003 R.Hordoir & Christian Ethe, Irregular |
---|
42 | ; grid on meridian axis |
---|
43 | ; |
---|
44 | ;- |
---|
45 | ;------------------------------------------------------------ |
---|
46 | ;------------------------------------------------------------ |
---|
47 | ;------------------------------------------------------------ |
---|
48 | PRO interp, znumdta, zdata_name,zsign, mask, t, glam, gphi, phi_field, zresul |
---|
49 | @common_interp |
---|
50 | tempsun = systime(1) ; pour key_performance |
---|
51 | ; |
---|
52 | ;------------------------------------------------- |
---|
53 | ; BILINEAR OR BICUBIC INTERPOLATION OF A FIELD * |
---|
54 | ;------------------------------------------------- |
---|
55 | ; |
---|
56 | ; 0. Extracting data at timestep t |
---|
57 | ; ================================ |
---|
58 | ; |
---|
59 | zzresul = fltarr(jpioce, jpjoce, jpkatm) |
---|
60 | zzresul2 = fltarr(jpioce, jpjoce, jpkoce) |
---|
61 | |
---|
62 | FOR jk = 0, jpkatm-1 DO BEGIN |
---|
63 | |
---|
64 | |
---|
65 | IF keyword_set(ndim) THEN $ |
---|
66 | ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1, 1], offset = [0, 0, jk, t] ELSE BEGIN |
---|
67 | IF keyword_set(key_2d) THEN $ |
---|
68 | ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ |
---|
69 | ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] |
---|
70 | ENDELSE |
---|
71 | |
---|
72 | zmask = fltarr(jpiatm, jpjatm) |
---|
73 | zmask = mask(*, *, jk) |
---|
74 | printf, 40, '' |
---|
75 | printf, 40, 'Data reading OK' |
---|
76 | IF keyword_set(key_vari) THEN BEGIN |
---|
77 | varid = ncdf_varid(znumdta, zdata_name) |
---|
78 | ncdf_attget, znumdta, varid, "scale_factor", scale_factor |
---|
79 | ncdf_attget, znumdta, varid, "add_offset", add_offset |
---|
80 | |
---|
81 | printf, 40, '' |
---|
82 | printf, 40, 'Scale factor and add_offset reading OK' |
---|
83 | ENDIF |
---|
84 | ; Scale Factor + Add_offset |
---|
85 | ; Modif Fred |
---|
86 | ; zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='cond_sed') |
---|
87 | zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='Br') |
---|
88 | |
---|
89 | ; Fin modif fred |
---|
90 | zdata = zdata*scale_factor + add_offset |
---|
91 | |
---|
92 | ; The data array and its mask are shifted zonally as defined in |
---|
93 | ; namelist by lon_shift |
---|
94 | |
---|
95 | zdata = shift(zdata, lon_shift, 0) |
---|
96 | zmask = shift(zmask, lon_shift, 0) |
---|
97 | |
---|
98 | ; zdata = double(zdata) |
---|
99 | IF keyword_set(nreverse) THEN BEGIN |
---|
100 | zdata = reverse(zdata, 2) |
---|
101 | printf, 40, 'Re-arranging data in latitudes : North-South --> South-North' |
---|
102 | ENDIF |
---|
103 | |
---|
104 | |
---|
105 | ; |
---|
106 | ; 0.5 Checking Coherence of winds on North Pole |
---|
107 | ; ============================================= |
---|
108 | ; |
---|
109 | |
---|
110 | |
---|
111 | |
---|
112 | IF keyword_set(nscal) AND keyword_set(north_pole) THEN BEGIN |
---|
113 | |
---|
114 | IF zdata_name EQ data_u_name then BEGIN |
---|
115 | |
---|
116 | |
---|
117 | zdata_u=zdata |
---|
118 | |
---|
119 | IF keyword_set(key_2d) THEN $ |
---|
120 | ncdf_varget, numdta_v, data_v_name, zdata_v, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ |
---|
121 | ncdf_varget, numdta_v, data_v_name, zdata_v, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] |
---|
122 | |
---|
123 | |
---|
124 | ENDIF |
---|
125 | |
---|
126 | IF zdata_name EQ data_v_name then BEGIN |
---|
127 | |
---|
128 | |
---|
129 | zdata_v=zdata |
---|
130 | |
---|
131 | IF keyword_set(key_2d) THEN $ |
---|
132 | ncdf_varget, numdta_u, data_u_name, zdata_u, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ |
---|
133 | ncdf_varget, numdta_u, data_u_name, zdata_u, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] |
---|
134 | |
---|
135 | |
---|
136 | ENDIF |
---|
137 | |
---|
138 | ; Check is made here |
---|
139 | |
---|
140 | global_vect=replicate(0.,jpiatm,jpjatm,2) |
---|
141 | |
---|
142 | global_vect(*,*,0)=zdata_u |
---|
143 | global_vect(*,*,1)=zdata_v |
---|
144 | |
---|
145 | northwind, global_vect,zdata_name,t |
---|
146 | |
---|
147 | zdata_u=global_vect(*,*,0) |
---|
148 | zdata_v=global_vect(*,*,1) |
---|
149 | |
---|
150 | ; End check and correction |
---|
151 | |
---|
152 | ; Puts back corrected value into initial variable |
---|
153 | |
---|
154 | |
---|
155 | IF zdata_name EQ data_u_name then BEGIN |
---|
156 | |
---|
157 | zdata=zdata_u |
---|
158 | |
---|
159 | ENDIF |
---|
160 | |
---|
161 | IF zdata_name EQ data_v_name then BEGIN |
---|
162 | |
---|
163 | zdata=zdata_v |
---|
164 | |
---|
165 | ENDIF |
---|
166 | |
---|
167 | ENDIF |
---|
168 | |
---|
169 | ; |
---|
170 | ; 1. Applying mask on data and plotting |
---|
171 | ; ===================================== |
---|
172 | ; |
---|
173 | val = 1.e20 |
---|
174 | IF keyword_set(nmiss) THEN BEGIN |
---|
175 | |
---|
176 | ind = where(abs(zdata) GE abs(missing_val)) |
---|
177 | IF ind[0] NE -1 then zmask(ind) = 0. |
---|
178 | |
---|
179 | ENDIF |
---|
180 | |
---|
181 | ind = where(zmask EQ 0.) |
---|
182 | |
---|
183 | IF ind[0] NE -1 then zdata(ind) = val |
---|
184 | |
---|
185 | ind = where(zmask NE 0.) |
---|
186 | IF ind[0] NE -1 THEN min = min(zdata(ind)) ELSE BEGIN |
---|
187 | min = missing_val |
---|
188 | ENDELSE |
---|
189 | IF ind[0] NE -1 THEN max = max(zdata(ind)) ELSE BEGIN |
---|
190 | max = missing_val |
---|
191 | zdata(*, *) = missing_val |
---|
192 | ENDELSE |
---|
193 | IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN |
---|
194 | IF keyword_set(key_ps) THEN openps, string(zdata_name, '_controle_iteration_', strtrim(t+1, 2), '.ps') |
---|
195 | draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 1, 3, $ |
---|
196 | TITLE = zdata_name+', iteration = '+strtrim(t+1,1), min, max |
---|
197 | ENDIF |
---|
198 | ; |
---|
199 | ; 2. Extrapolating coastal values on land points |
---|
200 | ; ============================================== |
---|
201 | ; |
---|
202 | eps = 0. |
---|
203 | n = 0 |
---|
204 | tempsun2 = systime(1) ; pour key_performance |
---|
205 | ; |
---|
206 | ; ... adding a border to overcome shifting problems |
---|
207 | ; |
---|
208 | zero = fltarr(jpiatm+2, jpjatm+2) |
---|
209 | zero[1:jpiatm, 1:jpjatm] = zmask |
---|
210 | zmask = zero |
---|
211 | zero = fltarr(jpiatm+2, jpjatm+2) |
---|
212 | zero[1:jpiatm, 1:jpjatm] = zdata |
---|
213 | zdata = zero |
---|
214 | |
---|
215 | |
---|
216 | ; |
---|
217 | ; ... filling |
---|
218 | ; |
---|
219 | |
---|
220 | if max(zmask) GT 0 then begin |
---|
221 | |
---|
222 | WHILE n LE nfil AND eps LE 1.e-5 DO BEGIN |
---|
223 | |
---|
224 | extrap, zdata, zmask, n, val |
---|
225 | eps = min(abs(zdata-replicate(val, jpiatm, jpjatm))) |
---|
226 | n = n+1 |
---|
227 | ENDWHILE |
---|
228 | |
---|
229 | endif |
---|
230 | ; |
---|
231 | ; ... recovering values on initial domain |
---|
232 | ; |
---|
233 | zdata = zdata[1:jpiatm, 1:jpjatm] |
---|
234 | zmask = zmask[1:jpiatm, 1:jpjatm] |
---|
235 | ; |
---|
236 | if keyword_set(key_performance) THEN print, 'temps extrap', systime(1)-tempsun2 |
---|
237 | ; .. drawing |
---|
238 | IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN |
---|
239 | draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 2, 3, $ |
---|
240 | TITLE = 'extrapolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max |
---|
241 | ENDIF |
---|
242 | |
---|
243 | ind = where(zmask NE 0. AND zdata GE val) |
---|
244 | IF n_elements(ind) GT 1 THEN BEGIN |
---|
245 | print, 'We stop because of insufficient extrapolation on missing values' |
---|
246 | print, n_elements(ind), ' points non extrapolated' |
---|
247 | print, 'Increase nfill in naminterp.pro' |
---|
248 | stop |
---|
249 | ENDIF |
---|
250 | |
---|
251 | |
---|
252 | ; 3. Treatment of Boundary Conditions |
---|
253 | ; ======================================================== |
---|
254 | ; Treatment of boundary conditions require extra stripes |
---|
255 | ; Size of extra stripe is calculated after a given extension |
---|
256 | ; along x and y axis on the OUTPUT grid. iplus and jplus will be |
---|
257 | ; the number of extra points required for the extra stripes of the |
---|
258 | ; INPUT grid. |
---|
259 | |
---|
260 | |
---|
261 | iplus=fix((ext_x*pres_x)/mean(dx)) |
---|
262 | jplus=fix((ext_y*pres_y)/mean(abs(phi_field(1:jpjatm/2)-shift(phi_field(1:jpjatm/2),1)))) |
---|
263 | |
---|
264 | ; Initialisation of array for the northern boundary condition |
---|
265 | |
---|
266 | zdata_north=replicate(0.,jpiatm,jpjatm+jplus) |
---|
267 | zdata_north(*,0:jpjatm-1)=zdata |
---|
268 | |
---|
269 | ; Data of the input grid is duplicated for the northern boundary condition |
---|
270 | |
---|
271 | zdata_inv=shift(zdata_north(*,jpjatm-1-jplus+1:jpjatm-1),jpiatm/2,0) |
---|
272 | |
---|
273 | zdata_inv=reverse(zdata_inv,2) |
---|
274 | |
---|
275 | zdata_north(*,jpjatm:jpjatm+jplus-1)=zdata_inv |
---|
276 | |
---|
277 | ; Transition to East/West Boundary Conditions |
---|
278 | |
---|
279 | zdata=zdata_north |
---|
280 | |
---|
281 | ; East/West Boundary Condition, zdata_extend will be the array that |
---|
282 | ; contains both northern boundary condition and E/W boundary |
---|
283 | ; condition, on the input grid |
---|
284 | |
---|
285 | |
---|
286 | zdata_extend=replicate(0.,jpiatm+2*iplus,jpjatm+jplus) |
---|
287 | |
---|
288 | zdata_extend(iplus:jpiatm-1+iplus,*)=zdata |
---|
289 | |
---|
290 | zdata_extend(0:iplus-1,*)=zdata(jpiatm-1-iplus+1:jpiatm-1,*) |
---|
291 | |
---|
292 | zdata_extend(jpiatm+iplus:jpiatm+2*iplus-1,*)=zdata(0:iplus-1,*) |
---|
293 | |
---|
294 | ; |
---|
295 | ; 3.1 Interpolating data on output grid and plotting result |
---|
296 | ; ========================================================= |
---|
297 | ; The interpolation is made bellow. We use the index of the input grid |
---|
298 | ; to interpolate which makes we need to convert the latitudes and |
---|
299 | ; longitudes of the output grid its dimensions. ext_x*pres_x/dx is added to |
---|
300 | ; take into account the extra stripe at the left handside of the input |
---|
301 | ; grid |
---|
302 | |
---|
303 | print, 'Time Step=',t |
---|
304 | print, 'Input Level=',jk |
---|
305 | |
---|
306 | |
---|
307 | glam2 = (glam - west_lon+ext_x*pres_x )/dx |
---|
308 | gphi2 = replicate(0.,jpioce,jpjoce) |
---|
309 | |
---|
310 | |
---|
311 | |
---|
312 | for ji = 0, jpioce - 1 do begin |
---|
313 | for jj = 0, jpjoce - 1 do begin |
---|
314 | |
---|
315 | diff = - gphi(ji,jj) $ |
---|
316 | + phi_field |
---|
317 | |
---|
318 | diff = where( diff GE 0. ) |
---|
319 | |
---|
320 | |
---|
321 | j_above=min( diff ) |
---|
322 | |
---|
323 | if diff(0) eq -1 then begin |
---|
324 | ; print, 'Output grid goes above input grid latitude top value !!!' |
---|
325 | ; print, 'Output grid top value has been take as input grid maximum value' |
---|
326 | ; print, 'Please check results at the vicinity of North Pole' |
---|
327 | |
---|
328 | j_fine = jpjatm-1 |
---|
329 | gphi2(ji,jj) = j_fine |
---|
330 | endif else begin |
---|
331 | |
---|
332 | j_bellow = j_above -1 |
---|
333 | |
---|
334 | if j_bellow GE 0 then begin |
---|
335 | |
---|
336 | j_fine = j_bellow |
---|
337 | |
---|
338 | if j_above LT jpjatm then begin |
---|
339 | j_fine = j_bellow+( gphi(ji,jj)-phi_field(j_bellow) )/( phi_field(j_above)-phi_field(j_bellow) ) |
---|
340 | endif |
---|
341 | |
---|
342 | gphi2(ji,jj) = j_fine |
---|
343 | if j_fine GT j_above then stop |
---|
344 | |
---|
345 | endif else begin |
---|
346 | |
---|
347 | j_fine = j_above |
---|
348 | |
---|
349 | gphi2(ji,jj) = j_fine |
---|
350 | |
---|
351 | endelse |
---|
352 | |
---|
353 | |
---|
354 | endelse |
---|
355 | |
---|
356 | endfor |
---|
357 | endfor |
---|
358 | |
---|
359 | |
---|
360 | IF keyword_set(ninterp) THEN BEGIN |
---|
361 | |
---|
362 | zresul = interpolate(zdata_extend, glam2, gphi2, cubic = -.5) |
---|
363 | |
---|
364 | ENDIF ELSE BEGIN |
---|
365 | |
---|
366 | zresul = bilinear(zdata_extend , glam2, gphi2) |
---|
367 | |
---|
368 | ENDELSE |
---|
369 | |
---|
370 | IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN |
---|
371 | draw, zresul, glam, gphi, jpioce, jpjoce, 1, 3, 3, $ |
---|
372 | TITLE = 'interpolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max |
---|
373 | IF keyword_set(key_ps) THEN BEGIN |
---|
374 | closeps |
---|
375 | ENDIF ELSE BEGIN |
---|
376 | print, 'Press return to continue' |
---|
377 | zxc = strarr(1) |
---|
378 | read, zxc |
---|
379 | ENDELSE |
---|
380 | ENDIF |
---|
381 | zzresul(*, *, jk) = temporary(zresul) |
---|
382 | ENDFOR |
---|
383 | zresul = zzresul |
---|
384 | |
---|
385 | ; |
---|
386 | ; 4. Interpolating data on the vertical axis if the outdept not equal datadept |
---|
387 | ; =========================================================================== |
---|
388 | |
---|
389 | |
---|
390 | IF keyword_set(ndim) THEN BEGIN |
---|
391 | |
---|
392 | nlev_flag = WHERE ( datadept(0:jpkatm-1) ne outdept(0:jpkoce-1) ) |
---|
393 | IF nlev_flag(0) NE -1 THEN BEGIN |
---|
394 | FOR jk = 0, jpkoce-1 DO BEGIN |
---|
395 | |
---|
396 | jn = where(abs(datadept) GE abs(outdept(jk)) $ |
---|
397 | AND abs(shift(datadept, 1)) LT abs(outdept(jk))) |
---|
398 | |
---|
399 | |
---|
400 | IF min(jn) GE 0 THEN BEGIN |
---|
401 | |
---|
402 | ikt=jn-1 |
---|
403 | ikb=jn |
---|
404 | |
---|
405 | IF jn NE 0 then begin |
---|
406 | |
---|
407 | zt1=abs(outdept(jk))-abs(datadept(ikt)) |
---|
408 | zt2=abs(datadept(ikb))-abs(outdept(jk)) |
---|
409 | ztx1=abs(datadept(ikb))-abs(datadept(ikt)) |
---|
410 | za1=replicate(zt1/ztx1, jpioce, jpjoce) |
---|
411 | zb1=replicate(zt2/ztx1, jpioce, jpjoce) |
---|
412 | zzresul2(*, *, jk) = temporary(za1*zzresul(*, *, ikb)+zb1*zzresul(*, *, ikt)) |
---|
413 | |
---|
414 | ENDIF ELSE BEGIN |
---|
415 | |
---|
416 | zzresul2(*, *, jk) = zzresul(*, *, ikb) |
---|
417 | |
---|
418 | ENDELSE |
---|
419 | |
---|
420 | ENDIF ELSE BEGIN |
---|
421 | |
---|
422 | zzresul2(*, *, jk) = temporary(zzresul2(*, *, jk-1)) |
---|
423 | |
---|
424 | ENDELSE |
---|
425 | ENDFOR |
---|
426 | zresul = zzresul2 |
---|
427 | ENDIF ELSE BEGIN |
---|
428 | print, ' ' |
---|
429 | print, 'No need to perform vertical interpolation' |
---|
430 | print, ' ' |
---|
431 | ENDELSE |
---|
432 | ENDIF |
---|
433 | |
---|
434 | |
---|
435 | if keyword_set(key_performance) THEN print, 'temps interp', systime(1)-tempsun |
---|
436 | |
---|
437 | return |
---|
438 | END |
---|