source: trunk/interp.pro @ 6

Last change on this file since 6 was 2, checked in by pinsard, 18 years ago

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