source: trunk/interpolation.pro @ 2

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

initial import from /usr/work/fvi/OPA/geomag/

File size: 12.4 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME: interpolation.pro
6;
7; PURPOSE: Interpolates data from a regular origin grid to any grid
8;
9; CATEGORY: Program
10;
11; CALLING SEQUENCE: interpolation
12;
13; INPUTS: Complete the INIT_PATH and NAMINTERP files     
14;
15; KEYWORD PARAMETERS: None
16;
17; OUTPUTS: NetCDF file(s) containing the interpolated field
18;
19; COMMON BLOCKS:
20;       common_interp.pro
21;
22; SIDE EFFECTS:
23;
24; RESTRICTIONS:
25;
26; EXAMPLE:
27;
28; MODIFICATION HISTORY: 11/99 A. Jouzeau
29;                       08/00 A. Jouzeau
30;                       01/02 R. Hordoir
31;                       03/2003 R. Hordoir & C. Ethe
32;
33;-
34;------------------------------------------------------------
35;------------------------------------------------------------
36;------------------------------------------------------------
37PRO interpolation
38@common_interp
39   tempsun = systime(1) ; pour key_performance
40print,'***********************************************'
41print,'********** OPA INTERPOLATION PACKAGE **********'
42print,'***********************************************'
43;
44unit = 40
45openw, unit, 'interp_output'
46;
47printf, 40, '*******************************************'
48printf, 40, '*         INTERPOLATION PROGRAM           *'
49printf, 40, '*******************************************'
50printf, 40, ''
51printf, 40, '******************************************'
52printf, 40, '* USE :                                  *'
53printf, 40, '* --> Input  grid : Regular              *'
54printf, 40, '* --> Output grid : Regular or Irregular *'
55printf, 40, '******************************************'
56;
57;--------------------------------------------------
58; INITIALISATIONS                                 *
59;--------------------------------------------------
60;
61printf, 40, ''
62printf, 40, 'I- Interpolation initialisations'
63printf, 40, '================================'
64;
65@naminterp
66;
67printf, 40, ''
68printf, 40, 'II- Paths initialisations'
69printf, 40, '========================='
70;
71@init_path
72
73;
74;---------------------------------------------------
75; READING INPUT MASK & OUTPUT COORDINATES          *
76;---------------------------------------------------
77;
78printf, 40, ''
79printf, 40, 'III- Opening data files and reading input mask and output coordinates'
80printf, 40, '====================================================================='
81;
82netcdf_input, mask
83;
84; Coordinates of input grid (for plotting purposes)
85;
86lonatm = fltarr(jpiatm)
87latatm = fltarr(jpjatm)
88FOR i = 0, jpiatm-1 DO BEGIN
89   lonatm(i) = west_lon + i*dx
90ENDFOR
91FOR j = 0, jpjatm-1 DO BEGIN
92   latatm(j) = phi_input(j)
93ENDFOR
94;
95;---------------------------------------------------
96; DATA MASK PREPROCESSING                          *
97;---------------------------------------------------
98;
99printf, 40, ''
100printf, 40, 'IV- Data mask preprocessing'
101printf, 40, '==========================='
102;
103IF keyword_set(negmask) THEN BEGIN
104   mask = abs(temporary(mask))
105   printf, 40, 'Taking the abs of mask (-1 --> +1)'
106ENDIF
107IF keyword_set(nreverse) THEN BEGIN
108   mask = reverse(mask, 2)
109   printf, 40, 'Re-arranging mask in latitudes : North-South --> South-North'
110ENDIF
111ind = where(mask NE 0 AND mask NE 1)
112IF n_elements(ind) NE 1 THEN BEGIN
113   printf, 40, ''
114   printf, 40, '************************************************************'
115   printf, 40, '* ERROR: data mask file must contain ONLY 0 and 1 values!! *'
116   printf, 40, '************************************************************'
117   printf, 40, ''
118   printf, 40, 'program stops'
119   stop
120ENDIF
121IF keyword_set(invmask) THEN BEGIN
122   mask = inv_mask(mask)
123   printf, 40, 'Giving value = 1 to mask on ocean points '
124ENDIF
125IF keyword_set(extend)  THEN BEGIN
126   mask = preproc_mask(mask)
127   printf, 40, 'Extending mask over one point on ocean points '
128ENDIF
129
130
131; mask = double(mask)
132
133;
134;---------------------------------------------------
135; INTERPOLATION                                    *
136;---------------------------------------------------
137;
138printf, 40, ''
139printf, 40, 'V- Interpolation'
140printf, 40, '================'
141;
142; 0. Creating output file and storing attributes
143; ==============================================
144;
145IF keyword_set(nscal) THEN BEGIN
146   printf, 40, ''
147   printf, 40, ' V-1- Writing attributes of interpolated field to files'
148   printf, 40, ' ------------------------------------------------------'
149;
150; ... x-component of interpolated field
151   printf, 40, ''
152   printf, 40, 'Writing attributes of x-component of interpolated field'
153   ncdf_varget, nummsh, mask_u_name, umask, count = [jpioce, jpjoce, jpkoce, 1]
154   netcdf_output, mask_u_name, umask, data_u_name, lon_u, lat_u, outputfile_x, $
155    title_x, long_x_name, POINT = "u", $
156    numout_u, varid_u, numdta_u
157; ... y-component of interpolated field
158   printf, 40, ''
159   printf, 40, 'Writing attributes of y-component of interpolated field'
160   ncdf_varget, nummsh, mask_v_name, vmask, count = [jpioce, jpjoce, jpkoce, 1]
161   netcdf_output, mask_v_name, vmask, data_v_name, lon_v, lat_v, outputfile_y, $
162    title_y, long_y_name, POINT = "v", $
163    numout_v, varid_v, numdta_v
164ENDIF ELSE BEGIN
165   printf, 40, ''
166   printf, 40, ' V-1- Writing attributes of interpolated field to file'
167   printf, 40, ' -----------------------------------------------------'
168;
169   ncdf_varget, nummsh, mask_t_name, tmask, count = [jpioce, jpjoce, jpkoce, 1]
170   netcdf_output, mask_t_name, tmask, data_name, lon_t, lat_t, outputfile, $
171    title, long_name, POINT = "t", $
172    numout, varid, numdta
173ENDELSE
174
175
176;
177; 1. Interpolation
178; ================
179;
180printf, 40, ''
181printf, 40, ' V-2- Interpolating, correcting and writing field(s)'
182printf, 40, ' ---------------------------------------------------'
183printf, 40, ''
184printf, 40, ' **** ENTERING TEMPORAL LOOP ****'
185;
186;  BEGINNING OF TIME LOOP
187;
188IF keyword_set(nscal) THEN BEGIN
189   mi_u = 1.e20
190   ma_u = -1.e20
191   mi_v = 1.e20
192   ma_v = -1.e20
193ENDIF ELSE BEGIN
194   mi = 1.e20
195   ma = -1.e20
196ENDELSE
197
198
199FOR t = nit000-1, nitend-1 DO BEGIN
200   tempsun4 = systime(1)        ; pour key_performance
201   printf, 40, ''
202   printf, 40, '  ***>> TIME ITERATION NUMBER :', t+1
203   IF keyword_set(ndim) THEN BEGIN
204      interp, numdta, data_name,1, mask, t, lon_t, lat_t, phi_input, zresul
205   ENDIF ELSE BEGIN
206      IF keyword_set(nscal) THEN BEGIN
207         printf, 40, 'A vectorial field has been chosen for interpolation'
208         IF keyword_set(north_pole) THEN BEGIN
209         printf, 40, 'Vectorial Fied will be Tested for Coherence nearby North Pole'
210         ENDIF
211         printf, 40, ''
212         printf, 40, 'interpolating ', data_u_name, ' at u-point and v-point '
213         interp, numdta_u, data_u_name,1, mask, t, lon_u, lat_u, phi_input_u, zresul_uu
214         interp, numdta_u, data_u_name,1, mask, t, lon_v, lat_v, phi_input_u, zresul_uv
215         printf, 40, 'interpolating ', data_v_name, ' at v-point and u-point '
216         interp, numdta_v, data_v_name,1 ,mask, t, lon_v, lat_v, phi_input_v, zresul_vv
217         interp, numdta_v, data_v_name,1,mask, t, lon_u, lat_u,  phi_input_v, zresul_vu
218      ENDIF ELSE BEGIN
219         printf, 40, ''
220         printf, 40, 'interpolating ', data_name, ' at t-point'
221         interp, numdta, data_name, 1,mask, t, lon_t, lat_t, phi_input, zresul
222      ENDELSE
223   ENDELSE
224
225
226
227;
228; 2. Angles correction for 2D fields                 
229; ==================================
230;
231   IF keyword_set(nscal) THEN BEGIN
232      printf, 40, ''
233      printf, 40, ' angles correction'
234;
235      tempsun2 = systime(1)      ; pour key_performance
236      correc_angle, zresul_uu, zresul_uv, zresul_vv, zresul_vu, zresul_u, zresul_v
237      if keyword_set(key_performance) THEN print, 'temps angles', systime(1)-tempsun2
238      IF t MOD ndraw EQ 0 AND keyword_set(keydraw) THEN BEGIN
239         ncdf_varget, numdta_u, data_u_name, zdata_u, $
240          count = [jpiatm, jpjatm, 1], offset = [0, 0, t]
241         ncdf_varget, numdta_v, data_v_name, zdata_v, $
242          count = [jpiatm, jpjatm, 1], offset = [0, 0, t]
243         zdata_u = zdata_u*scale_factor + add_offset
244         zdata_v = zdata_v*scale_factor + add_offset
245         IF keyword_set(nreverse) THEN zdata_u = reverse(zdata_u, 2)
246         IF keyword_set(nreverse) THEN zdata_v = reverse(zdata_v, 2)
247         zu = zdata_u
248         zv = zdata_v
249         IF keyword_set(nmiss) THEN BEGIN
250            ind = where(abs(zu) GE missing_val)
251            zu(ind) = 0.
252            ind = where(abs(zv) GE missing_val)
253            zv(ind) = 0.
254         ENDIF
255         ind = where(mask NE 0.)
256         minu = min(zu(ind))
257         maxu = max(zu(ind))
258         minv = min(zv(ind))
259         maxv = max(zv(ind))
260         IF keyword_set(key_ps) THEN $
261          openps, string('control_corrected_winds_iteration_', strtrim(t+1, 2), '.ps')
262         draw, zdata_u, lonatm, latatm, jpiatm, jpjatm, 2, 1, 2, $
263          TITLE = data_u_name+' data, iteration = '+strtrim(t+1,1), minu, maxu, /land
264         draw, zdata_v, lonatm, latatm, jpiatm, jpjatm, 2, 2, 2, $
265          TITLE = data_v_name+' data, iteration = '+strtrim(t+1,1), minv, maxv, /land
266         draw, zresul_u, lon_f, lat_f, jpioce, jpjoce, 2, 3, 2, $
267          TITLE = data_u_name+' interpolated, iteration = '+strtrim(t+1,1), minu, maxu, /land
268         draw, zresul_v, lon_f, lat_f, jpioce, jpjoce, 2, 4, 2, $
269          TITLE = data_v_name+' interpolated, iteration = '+strtrim(t+1,1), minv, maxv, /land
270         IF keyword_set(key_ps) THEN BEGIN
271            closeps
272         ENDIF ELSE BEGIN
273            print, 'Press return to continue'
274            zxc = strarr(1)
275            read, zxc
276         ENDELSE
277      ENDIF
278   ENDIF
279 
280;
281; 3. NetCDF output to file of interpolated field
282; ==============================================
283;
284    tempsun3 = systime(1) ; pour key_performance
285  IF keyword_set(nscal) THEN BEGIN
286      printf, 40, ''
287      printf, 40, 'Storing components of interpolated field'
288     
289
290      ncdf_varput, numout_u, varid_u, zresul_u, offset = [0, 0, 0, t-nit000+1]
291      ind = where(abs(zresul_u) LT 1.e10)
292      IF ind[0] NE -1 THEN BEGIN
293         IF min(zresul_u) LT mi_u THEN mi_u = min(zresul_u[ind])
294         IF max(zresul_u) GT ma_u THEN ma_u = max(zresul_u[ind])
295      ENDIF ELSE BEGIN
296         mi_u = missing_val
297         ma_u = missing_val
298      ENDELSE
299       
300      ncdf_varput, numout_v, varid_v, zresul_v, offset = [0, 0, 0, t-nit000+1]
301      ind = where(abs(zresul_v) LT 1.e10)
302      IF ind[0] NE -1 THEN BEGIN
303         IF min(zresul_v) LT mi_v THEN mi_v = min(zresul_v[ind])
304         IF max(zresul_v) GT ma_v THEN ma_v = max(zresul_v[ind])
305      ENDIF ELSE BEGIN
306         mi_v = missing_val
307         ma_v = missing_val
308      ENDELSE
309   ENDIF ELSE BEGIN
310      printf, 40, ''
311      printf, 40, 'Storing field'
312     
313      ncdf_varput, numout, varid, zresul, offset = [0, 0, 0, t-nit000+1]
314      ind = where(abs(zresul) LT 1.e10)
315      IF min(zresul) LT mi THEN mi = min(zresul[ind])
316      IF max(zresul) GT ma THEN ma = max(zresul[ind])
317   ENDELSE
318      if keyword_set(key_performance) THEN print, 'temps ecriture variable', systime(1)-tempsun3
319      if keyword_set(key_performance) THEN printf, 40, 'temps ecriture variable', systime(1)-tempsun3
320      if keyword_set(key_performance) THEN print, 'temps boucle', systime(1)-tempsun4
321      if keyword_set(key_performance) THEN printf, 40, 'temps boucle', systime(1)-tempsun4
322ENDFOR
323;
324;  END OF TIME LOOP
325;
326printf, 40, ''
327printf, 40, '**** END OF TEMPORAL LOOP ****'
328;
329printf, 40, ''
330printf, 40, 'Closing files'
331;
332; storing min and max in attributes
333;
334IF keyword_set(nscal) THEN BEGIN
335   ncdf_control, numout_u, /redef
336   ncdf_control, numout_v, /redef
337   ncdf_attput, numout_u, varid_u, "valid_min", mi_u
338   ncdf_attput, numout_u, varid_u, "valid_max", ma_u
339   ncdf_attput, numout_v, varid_v, "valid_min", mi_v
340   ncdf_attput, numout_v, varid_v, "valid_max", ma_v
341ENDIF ELSE BEGIN
342   ncdf_control, numout, /redef
343   ncdf_attput, numout, varid, "valid_min", mi
344   ncdf_attput, numout, varid, "valid_max", ma
345ENDELSE
346IF keyword_set(key_mask) THEN ncdf_close, nummsk
347IF keyword_set(nscal) THEN BEGIN
348   ncdf_close, numdta_u
349   ncdf_close, numdta_v
350   ncdf_close, numout_u
351   ncdf_close, numout_v
352ENDIF ELSE BEGIN
353   ncdf_close, numout
354   ncdf_close, numdta
355ENDELSE
356ncdf_close, nummsh
357;
358printf, 40, ''
359printf, 40, '*******************************************'
360printf, 40, '*      END OF INTERPOLATION PROGRAM       *'
361printf, 40, '*******************************************'
362;
363close, 40
364   if keyword_set(key_performance) THEN print, 'temps total', systime(1)-tempsun
365   if keyword_set(key_performance) THEN printf,40, 'temps total', systime(1)-tempsun
366;
367return
368END
Note: See TracBrowser for help on using the repository browser.