;+ ; NAME: interp.pro ; ; PURPOSE: Performs a bilinear or bicubic interpolation on ; data from a regular origin grid to any grid ; CATEGORY: Subroutine ; ; CALLING SEQUENCE: interp, znumdta, zdata_name, mask, t, glam, gphi, zresul ; ; INPUTS: ; ; znumdta : ID of input data file ; zdata_name : input data name ; mask : data mask ; t : timestep ; glam : longitude of target grid ; gphi : latitude " " " ; ; KEYWORD PARAMETERS: None ; ; OUTPUTS: ; zresul : interpolated field ; ; COMMON BLOCKS: ; common_interp.pro ; ; SIDE EFFECTS: ; ; RESTRICTIONS: ; ; EXAMPLE: ; ; MODIFICATION HISTORY: 11/99 A. Jouzeau ; 04/00 A. Jouzeau, M. Imbard ; 08/00 A. Jouzeau ; 10/01 R. Hordoir : Added boundary conditions, ; Bug correction, North Pole Problem ; 02/02 R. Hordoir : Added North Pole Treatment ; for non coherent winds. ; 03/2003 R.Hordoir & Christian Ethe, Irregular ; grid on meridian axis ; ;- ;------------------------------------------------------------ ;------------------------------------------------------------ ;------------------------------------------------------------ PRO interp, znumdta, zdata_name,zsign, mask, t, glam, gphi, phi_field, zresul @common_interp tempsun = systime(1) ; pour key_performance ; ;------------------------------------------------- ; BILINEAR OR BICUBIC INTERPOLATION OF A FIELD * ;------------------------------------------------- ; ; 0. Extracting data at timestep t ; ================================ ; zzresul = fltarr(jpioce, jpjoce, jpkatm) zzresul2 = fltarr(jpioce, jpjoce, jpkoce) FOR jk = 0, jpkatm-1 DO BEGIN IF keyword_set(ndim) THEN $ ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1, 1], offset = [0, 0, jk, t] ELSE BEGIN IF keyword_set(key_2d) THEN $ ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ ncdf_varget, znumdta, zdata_name, zdata, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] ENDELSE zmask = fltarr(jpiatm, jpjatm) zmask = mask(*, *, jk) printf, 40, '' printf, 40, 'Data reading OK' IF keyword_set(key_vari) THEN BEGIN varid = ncdf_varid(znumdta, zdata_name) ncdf_attget, znumdta, varid, "scale_factor", scale_factor ncdf_attget, znumdta, varid, "add_offset", add_offset printf, 40, '' printf, 40, 'Scale factor and add_offset reading OK' ENDIF ; Scale Factor + Add_offset ; Modif Fred ; zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='cond_sed') zdata=ncdf_lec('/usr/work/sur/fvi/OPA/geomag/condmag.nc',var='Br') ; Fin modif fred zdata = zdata*scale_factor + add_offset ; The data array and its mask are shifted zonally as defined in ; namelist by lon_shift zdata = shift(zdata, lon_shift, 0) zmask = shift(zmask, lon_shift, 0) ; zdata = double(zdata) IF keyword_set(nreverse) THEN BEGIN zdata = reverse(zdata, 2) printf, 40, 'Re-arranging data in latitudes : North-South --> South-North' ENDIF ; ; 0.5 Checking Coherence of winds on North Pole ; ============================================= ; IF keyword_set(nscal) AND keyword_set(north_pole) THEN BEGIN IF zdata_name EQ data_u_name then BEGIN zdata_u=zdata IF keyword_set(key_2d) THEN $ ncdf_varget, numdta_v, data_v_name, zdata_v, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ ncdf_varget, numdta_v, data_v_name, zdata_v, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] ENDIF IF zdata_name EQ data_v_name then BEGIN zdata_v=zdata IF keyword_set(key_2d) THEN $ ncdf_varget, numdta_u, data_u_name, zdata_u, count = [jpiatm, jpjatm, 1], offset = [0, 0, t] ELSE $ ncdf_varget, numdta_u, data_u_name, zdata_u, count = [jpiatm, jpjatm,1, 1], offset = [0, 0, 0, t] ENDIF ; Check is made here global_vect=replicate(0.,jpiatm,jpjatm,2) global_vect(*,*,0)=zdata_u global_vect(*,*,1)=zdata_v northwind, global_vect,zdata_name,t zdata_u=global_vect(*,*,0) zdata_v=global_vect(*,*,1) ; End check and correction ; Puts back corrected value into initial variable IF zdata_name EQ data_u_name then BEGIN zdata=zdata_u ENDIF IF zdata_name EQ data_v_name then BEGIN zdata=zdata_v ENDIF ENDIF ; ; 1. Applying mask on data and plotting ; ===================================== ; val = 1.e20 IF keyword_set(nmiss) THEN BEGIN ind = where(abs(zdata) GE abs(missing_val)) IF ind[0] NE -1 then zmask(ind) = 0. ENDIF ind = where(zmask EQ 0.) IF ind[0] NE -1 then zdata(ind) = val ind = where(zmask NE 0.) IF ind[0] NE -1 THEN min = min(zdata(ind)) ELSE BEGIN min = missing_val ENDELSE IF ind[0] NE -1 THEN max = max(zdata(ind)) ELSE BEGIN max = missing_val zdata(*, *) = missing_val ENDELSE IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN IF keyword_set(key_ps) THEN openps, string(zdata_name, '_controle_iteration_', strtrim(t+1, 2), '.ps') draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 1, 3, $ TITLE = zdata_name+', iteration = '+strtrim(t+1,1), min, max ENDIF ; ; 2. Extrapolating coastal values on land points ; ============================================== ; eps = 0. n = 0 tempsun2 = systime(1) ; pour key_performance ; ; ... adding a border to overcome shifting problems ; zero = fltarr(jpiatm+2, jpjatm+2) zero[1:jpiatm, 1:jpjatm] = zmask zmask = zero zero = fltarr(jpiatm+2, jpjatm+2) zero[1:jpiatm, 1:jpjatm] = zdata zdata = zero ; ; ... filling ; if max(zmask) GT 0 then begin WHILE n LE nfil AND eps LE 1.e-5 DO BEGIN extrap, zdata, zmask, n, val eps = min(abs(zdata-replicate(val, jpiatm, jpjatm))) n = n+1 ENDWHILE endif ; ; ... recovering values on initial domain ; zdata = zdata[1:jpiatm, 1:jpjatm] zmask = zmask[1:jpiatm, 1:jpjatm] ; if keyword_set(key_performance) THEN print, 'temps extrap', systime(1)-tempsun2 ; .. drawing IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN draw, zdata, lonatm, latatm, jpiatm, jpjatm, 1, 2, 3, $ TITLE = 'extrapolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max ENDIF ind = where(zmask NE 0. AND zdata GE val) IF n_elements(ind) GT 1 THEN BEGIN print, 'We stop because of insufficient extrapolation on missing values' print, n_elements(ind), ' points non extrapolated' print, 'Increase nfill in naminterp.pro' stop ENDIF ; 3. Treatment of Boundary Conditions ; ======================================================== ; Treatment of boundary conditions require extra stripes ; Size of extra stripe is calculated after a given extension ; along x and y axis on the OUTPUT grid. iplus and jplus will be ; the number of extra points required for the extra stripes of the ; INPUT grid. iplus=fix((ext_x*pres_x)/mean(dx)) jplus=fix((ext_y*pres_y)/mean(abs(phi_field(1:jpjatm/2)-shift(phi_field(1:jpjatm/2),1)))) ; Initialisation of array for the northern boundary condition zdata_north=replicate(0.,jpiatm,jpjatm+jplus) zdata_north(*,0:jpjatm-1)=zdata ; Data of the input grid is duplicated for the northern boundary condition zdata_inv=shift(zdata_north(*,jpjatm-1-jplus+1:jpjatm-1),jpiatm/2,0) zdata_inv=reverse(zdata_inv,2) zdata_north(*,jpjatm:jpjatm+jplus-1)=zdata_inv ; Transition to East/West Boundary Conditions zdata=zdata_north ; East/West Boundary Condition, zdata_extend will be the array that ; contains both northern boundary condition and E/W boundary ; condition, on the input grid zdata_extend=replicate(0.,jpiatm+2*iplus,jpjatm+jplus) zdata_extend(iplus:jpiatm-1+iplus,*)=zdata zdata_extend(0:iplus-1,*)=zdata(jpiatm-1-iplus+1:jpiatm-1,*) zdata_extend(jpiatm+iplus:jpiatm+2*iplus-1,*)=zdata(0:iplus-1,*) ; ; 3.1 Interpolating data on output grid and plotting result ; ========================================================= ; The interpolation is made bellow. We use the index of the input grid ; to interpolate which makes we need to convert the latitudes and ; longitudes of the output grid its dimensions. ext_x*pres_x/dx is added to ; take into account the extra stripe at the left handside of the input ; grid print, 'Time Step=',t print, 'Input Level=',jk glam2 = (glam - west_lon+ext_x*pres_x )/dx gphi2 = replicate(0.,jpioce,jpjoce) for ji = 0, jpioce - 1 do begin for jj = 0, jpjoce - 1 do begin diff = - gphi(ji,jj) $ + phi_field diff = where( diff GE 0. ) j_above=min( diff ) if diff(0) eq -1 then begin ; print, 'Output grid goes above input grid latitude top value !!!' ; print, 'Output grid top value has been take as input grid maximum value' ; print, 'Please check results at the vicinity of North Pole' j_fine = jpjatm-1 gphi2(ji,jj) = j_fine endif else begin j_bellow = j_above -1 if j_bellow GE 0 then begin j_fine = j_bellow if j_above LT jpjatm then begin j_fine = j_bellow+( gphi(ji,jj)-phi_field(j_bellow) )/( phi_field(j_above)-phi_field(j_bellow) ) endif gphi2(ji,jj) = j_fine if j_fine GT j_above then stop endif else begin j_fine = j_above gphi2(ji,jj) = j_fine endelse endelse endfor endfor IF keyword_set(ninterp) THEN BEGIN zresul = interpolate(zdata_extend, glam2, gphi2, cubic = -.5) ENDIF ELSE BEGIN zresul = bilinear(zdata_extend , glam2, gphi2) ENDELSE IF t MOD ndraw EQ 0 AND keyword_set(keydraw) AND jk MOD nlevel EQ 0 THEN BEGIN draw, zresul, glam, gphi, jpioce, jpjoce, 1, 3, 3, $ TITLE = 'interpolated '+zdata_name+', iteration = '+strtrim(t+1,1), min, max IF keyword_set(key_ps) THEN BEGIN closeps ENDIF ELSE BEGIN print, 'Press return to continue' zxc = strarr(1) read, zxc ENDELSE ENDIF zzresul(*, *, jk) = temporary(zresul) ENDFOR zresul = zzresul ; ; 4. Interpolating data on the vertical axis if the outdept not equal datadept ; =========================================================================== IF keyword_set(ndim) THEN BEGIN nlev_flag = WHERE ( datadept(0:jpkatm-1) ne outdept(0:jpkoce-1) ) IF nlev_flag(0) NE -1 THEN BEGIN FOR jk = 0, jpkoce-1 DO BEGIN jn = where(abs(datadept) GE abs(outdept(jk)) $ AND abs(shift(datadept, 1)) LT abs(outdept(jk))) IF min(jn) GE 0 THEN BEGIN ikt=jn-1 ikb=jn IF jn NE 0 then begin zt1=abs(outdept(jk))-abs(datadept(ikt)) zt2=abs(datadept(ikb))-abs(outdept(jk)) ztx1=abs(datadept(ikb))-abs(datadept(ikt)) za1=replicate(zt1/ztx1, jpioce, jpjoce) zb1=replicate(zt2/ztx1, jpioce, jpjoce) zzresul2(*, *, jk) = temporary(za1*zzresul(*, *, ikb)+zb1*zzresul(*, *, ikt)) ENDIF ELSE BEGIN zzresul2(*, *, jk) = zzresul(*, *, ikb) ENDELSE ENDIF ELSE BEGIN zzresul2(*, *, jk) = temporary(zzresul2(*, *, jk-1)) ENDELSE ENDFOR zresul = zzresul2 ENDIF ELSE BEGIN print, ' ' print, 'No need to perform vertical interpolation' print, ' ' ENDELSE ENDIF if keyword_set(key_performance) THEN print, 'temps interp', systime(1)-tempsun return END