source: trunk/interp.pro @ 48

Last change on this file since 48 was 48, checked in by pinsard, 10 years ago

fix thanks to coding rules

File size: 11.9 KB
Line 
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;------------------------------------------------------------
48PRO 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;
59zzresul = fltarr(jpioce, jpjoce, jpkatm)
60zzresul2 = fltarr(jpioce, jpjoce, jpkoce)
61
62FOR 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)
382ENDFOR
383zresul = zzresul
384
385;
386; 4. Interpolating data on the vertical axis if the outdept not equal datadept
387; ===========================================================================
388
389
390IF 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
432ENDIF
433
434
435   if keyword_set(key_performance) THEN print, 'temps interp', systime(1)-tempsun
436
437return
438END
Note: See TracBrowser for help on using the repository browser.