source: trunk/interpolation.pro

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

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