;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ ;+ ; NAME: interpolation.pro ; ; PURPOSE: Interpolates data from a regular origin grid to any grid ; ; CATEGORY: Program ; ; CALLING SEQUENCE: interpolation ; ; INPUTS: Complete the INIT_PATH and NAMINTERP files ; ; KEYWORD PARAMETERS: None ; ; OUTPUTS: NetCDF file(s) containing the interpolated field ; ; COMMON BLOCKS: ; common_interp.pro ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: 11/99 A. Jouzeau ; 08/00 A. Jouzeau ; 01/02 R. Hordoir ; 03/2003 R. Hordoir & C. Ethe ; ;- ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ PRO interpolation @common_interp tempsun = systime(1) ; pour key_performance print,'***********************************************' print,'********** OPA INTERPOLATION PACKAGE **********' print,'***********************************************' ; unit = 40 openw, unit, 'interp_output' ; printf, 40, '*******************************************' printf, 40, '* INTERPOLATION PROGRAM *' printf, 40, '*******************************************' printf, 40, '' printf, 40, '******************************************' printf, 40, '* USE : *' printf, 40, '* --> Input grid : Regular *' printf, 40, '* --> Output grid : Regular or Irregular *' printf, 40, '******************************************' ; ;-------------------------------------------------- ; INITIALISATIONS * ;-------------------------------------------------- ; printf, 40, '' printf, 40, 'I- Interpolation initialisations' printf, 40, '================================' ; @naminterp ; printf, 40, '' printf, 40, 'II- Paths initialisations' printf, 40, '=========================' ; @init_path ; ;--------------------------------------------------- ; READING INPUT MASK & OUTPUT COORDINATES * ;--------------------------------------------------- ; printf, 40, '' printf, 40, 'III- Opening data files and reading input mask and output coordinates' printf, 40, '=====================================================================' ; netcdf_input, mask ; ; Coordinates of input grid (for plotting purposes) ; lonatm = fltarr(jpiatm) latatm = fltarr(jpjatm) FOR i = 0, jpiatm-1 DO BEGIN lonatm(i) = west_lon + i*dx ENDFOR FOR j = 0, jpjatm-1 DO BEGIN latatm(j) = phi_input(j) ENDFOR ; ;--------------------------------------------------- ; DATA MASK PREPROCESSING * ;--------------------------------------------------- ; printf, 40, '' printf, 40, 'IV- Data mask preprocessing' printf, 40, '===========================' ; IF keyword_set(negmask) THEN BEGIN mask = abs(temporary(mask)) printf, 40, 'Taking the abs of mask (-1 --> +1)' ENDIF IF keyword_set(nreverse) THEN BEGIN mask = reverse(mask, 2) printf, 40, 'Re-arranging mask in latitudes : North-South --> South-North' ENDIF ind = where(mask NE 0 AND mask NE 1) IF n_elements(ind) NE 1 THEN BEGIN printf, 40, '' printf, 40, '************************************************************' printf, 40, '* ERROR: data mask file must contain ONLY 0 and 1 values!! *' printf, 40, '************************************************************' printf, 40, '' printf, 40, 'program stops' stop ENDIF IF keyword_set(invmask) THEN BEGIN mask = inv_mask(mask) printf, 40, 'Giving value = 1 to mask on ocean points ' ENDIF IF keyword_set(extend) THEN BEGIN mask = preproc_mask(mask) printf, 40, 'Extending mask over one point on ocean points ' ENDIF ; mask = double(mask) ; ;--------------------------------------------------- ; INTERPOLATION * ;--------------------------------------------------- ; printf, 40, '' printf, 40, 'V- Interpolation' printf, 40, '================' ; ; 0. Creating output file and storing attributes ; ============================================== ; IF keyword_set(nscal) THEN BEGIN printf, 40, '' printf, 40, ' V-1- Writing attributes of interpolated field to files' printf, 40, ' ------------------------------------------------------' ; ; ... x-component of interpolated field printf, 40, '' printf, 40, 'Writing attributes of x-component of interpolated field' ncdf_varget, nummsh, mask_u_name, umask, count = [jpioce, jpjoce, jpkoce, 1] netcdf_output, mask_u_name, umask, data_u_name, lon_u, lat_u, outputfile_x, $ title_x, long_x_name, POINT = "u", $ numout_u, varid_u, numdta_u ; ... y-component of interpolated field printf, 40, '' printf, 40, 'Writing attributes of y-component of interpolated field' ncdf_varget, nummsh, mask_v_name, vmask, count = [jpioce, jpjoce, jpkoce, 1] netcdf_output, mask_v_name, vmask, data_v_name, lon_v, lat_v, outputfile_y, $ title_y, long_y_name, POINT = "v", $ numout_v, varid_v, numdta_v ENDIF ELSE BEGIN printf, 40, '' printf, 40, ' V-1- Writing attributes of interpolated field to file' printf, 40, ' -----------------------------------------------------' ; ncdf_varget, nummsh, mask_t_name, tmask, count = [jpioce, jpjoce, jpkoce, 1] netcdf_output, mask_t_name, tmask, data_name, lon_t, lat_t, outputfile, $ title, long_name, POINT = "t", $ numout, varid, numdta ENDELSE ; ; 1. Interpolation ; ================ ; printf, 40, '' printf, 40, ' V-2- Interpolating, correcting and writing field(s)' printf, 40, ' ---------------------------------------------------' printf, 40, '' printf, 40, ' **** ENTERING TEMPORAL LOOP ****' ; ; BEGINNING OF TIME LOOP ; IF keyword_set(nscal) THEN BEGIN mi_u = 1.e20 ma_u = -1.e20 mi_v = 1.e20 ma_v = -1.e20 ENDIF ELSE BEGIN mi = 1.e20 ma = -1.e20 ENDELSE FOR t = nit000-1, nitend-1 DO BEGIN tempsun4 = systime(1) ; pour key_performance printf, 40, '' printf, 40, ' ***>> TIME ITERATION NUMBER :', t+1 IF keyword_set(ndim) THEN BEGIN interp, numdta, data_name,1, mask, t, lon_t, lat_t, phi_input, zresul ENDIF ELSE BEGIN IF keyword_set(nscal) THEN BEGIN printf, 40, 'A vectorial field has been chosen for interpolation' IF keyword_set(north_pole) THEN BEGIN printf, 40, 'Vectorial Fied will be Tested for Coherence nearby North Pole' ENDIF printf, 40, '' printf, 40, 'interpolating ', data_u_name, ' at u-point and v-point ' interp, numdta_u, data_u_name,1, mask, t, lon_u, lat_u, phi_input_u, zresul_uu interp, numdta_u, data_u_name,1, mask, t, lon_v, lat_v, phi_input_u, zresul_uv printf, 40, 'interpolating ', data_v_name, ' at v-point and u-point ' interp, numdta_v, data_v_name,1 ,mask, t, lon_v, lat_v, phi_input_v, zresul_vv interp, numdta_v, data_v_name,1,mask, t, lon_u, lat_u, phi_input_v, zresul_vu ENDIF ELSE BEGIN printf, 40, '' printf, 40, 'interpolating ', data_name, ' at t-point' interp, numdta, data_name, 1,mask, t, lon_t, lat_t, phi_input, zresul ENDELSE ENDELSE ; ; 2. Angles correction for 2D fields ; ================================== ; IF keyword_set(nscal) THEN BEGIN printf, 40, '' printf, 40, ' angles correction' ; tempsun2 = systime(1) ; pour key_performance correc_angle, zresul_uu, zresul_uv, zresul_vv, zresul_vu, zresul_u, zresul_v if keyword_set(key_performance) THEN print, 'temps angles', systime(1)-tempsun2 IF t MOD ndraw EQ 0 AND keyword_set(keydraw) THEN BEGIN ncdf_varget, numdta_u, data_u_name, zdata_u, $ count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ncdf_varget, numdta_v, data_v_name, zdata_v, $ count = [jpiatm, jpjatm, 1], offset = [0, 0, t] zdata_u = zdata_u*scale_factor + add_offset zdata_v = zdata_v*scale_factor + add_offset IF keyword_set(nreverse) THEN zdata_u = reverse(zdata_u, 2) IF keyword_set(nreverse) THEN zdata_v = reverse(zdata_v, 2) zu = zdata_u zv = zdata_v IF keyword_set(nmiss) THEN BEGIN ind = where(abs(zu) GE missing_val) zu(ind) = 0. ind = where(abs(zv) GE missing_val) zv(ind) = 0. ENDIF ind = where(mask NE 0.) minu = min(zu(ind)) maxu = max(zu(ind)) minv = min(zv(ind)) maxv = max(zv(ind)) IF keyword_set(key_ps) THEN $ openps, string('control_corrected_winds_iteration_', strtrim(t+1, 2), '.ps') draw, zdata_u, lonatm, latatm, jpiatm, jpjatm, 2, 1, 2, $ TITLE = data_u_name+' data, iteration = '+strtrim(t+1,1), minu, maxu, /land draw, zdata_v, lonatm, latatm, jpiatm, jpjatm, 2, 2, 2, $ TITLE = data_v_name+' data, iteration = '+strtrim(t+1,1), minv, maxv, /land draw, zresul_u, lon_f, lat_f, jpioce, jpjoce, 2, 3, 2, $ TITLE = data_u_name+' interpolated, iteration = '+strtrim(t+1,1), minu, maxu, /land draw, zresul_v, lon_f, lat_f, jpioce, jpjoce, 2, 4, 2, $ TITLE = data_v_name+' interpolated, iteration = '+strtrim(t+1,1), minv, maxv, /land IF keyword_set(key_ps) THEN BEGIN closeps ENDIF ELSE BEGIN print, 'Press return to continue' zxc = strarr(1) read, zxc ENDELSE ENDIF ENDIF ; ; 3. NetCDF output to file of interpolated field ; ============================================== ; tempsun3 = systime(1) ; pour key_performance IF keyword_set(nscal) THEN BEGIN printf, 40, '' printf, 40, 'Storing components of interpolated field' ncdf_varput, numout_u, varid_u, zresul_u, offset = [0, 0, 0, t-nit000+1] ind = where(abs(zresul_u) LT 1.e10) IF ind[0] NE -1 THEN BEGIN IF min(zresul_u) LT mi_u THEN mi_u = min(zresul_u[ind]) IF max(zresul_u) GT ma_u THEN ma_u = max(zresul_u[ind]) ENDIF ELSE BEGIN mi_u = missing_val ma_u = missing_val ENDELSE ncdf_varput, numout_v, varid_v, zresul_v, offset = [0, 0, 0, t-nit000+1] ind = where(abs(zresul_v) LT 1.e10) IF ind[0] NE -1 THEN BEGIN IF min(zresul_v) LT mi_v THEN mi_v = min(zresul_v[ind]) IF max(zresul_v) GT ma_v THEN ma_v = max(zresul_v[ind]) ENDIF ELSE BEGIN mi_v = missing_val ma_v = missing_val ENDELSE ENDIF ELSE BEGIN printf, 40, '' printf, 40, 'Storing field' ncdf_varput, numout, varid, zresul, offset = [0, 0, 0, t-nit000+1] ind = where(abs(zresul) LT 1.e10) IF min(zresul) LT mi THEN mi = min(zresul[ind]) IF max(zresul) GT ma THEN ma = max(zresul[ind]) ENDELSE if keyword_set(key_performance) THEN print, 'temps ecriture variable', systime(1)-tempsun3 if keyword_set(key_performance) THEN printf, 40, 'temps ecriture variable', systime(1)-tempsun3 if keyword_set(key_performance) THEN print, 'temps boucle', systime(1)-tempsun4 if keyword_set(key_performance) THEN printf, 40, 'temps boucle', systime(1)-tempsun4 ENDFOR ; ; END OF TIME LOOP ; printf, 40, '' printf, 40, '**** END OF TEMPORAL LOOP ****' ; printf, 40, '' printf, 40, 'Closing files' ; ; storing min and max in attributes ; IF keyword_set(nscal) THEN BEGIN ncdf_control, numout_u, /redef ncdf_control, numout_v, /redef ncdf_attput, numout_u, varid_u, "valid_min", mi_u ncdf_attput, numout_u, varid_u, "valid_max", ma_u ncdf_attput, numout_v, varid_v, "valid_min", mi_v ncdf_attput, numout_v, varid_v, "valid_max", ma_v ENDIF ELSE BEGIN ncdf_control, numout, /redef ncdf_attput, numout, varid, "valid_min", mi ncdf_attput, numout, varid, "valid_max", ma ENDELSE IF keyword_set(key_mask) THEN ncdf_close, nummsk IF keyword_set(nscal) THEN BEGIN ncdf_close, numdta_u ncdf_close, numdta_v ncdf_close, numout_u ncdf_close, numout_v ENDIF ELSE BEGIN ncdf_close, numout ncdf_close, numdta ENDELSE ncdf_close, nummsh ; printf, 40, '' printf, 40, '*******************************************' printf, 40, '* END OF INTERPOLATION PROGRAM *' printf, 40, '*******************************************' ; close, 40 if keyword_set(key_performance) THEN print, 'temps total', systime(1)-tempsun if keyword_set(key_performance) THEN printf,40, 'temps total', systime(1)-tempsun ; return END