PRO Read_reference_CFS, input_file, $ pressure, $ temperature, $ specific, $ ; specific humidity vapor_pressure, $ ; for refractivity Geopotential_height ncfileid = NCDF_OPEN(input_file) NCDF_VARGET, ncfileid, 'TMP_P0_L100_GLL0', temperature ; K, [720, 361, 37], float NCDF_VARGET, ncfileid, 'SPFH_P0_L100_GLL0', specific; kg/kg, [720, 361, 37], float NCDF_VARGET, ncfileid, 'HGT_P0_L100_GLL0', Geopotential_height ; gpm, [720, 361, 37], float NCDF_CLOSE, ncfileid temperature = double(temperature) specific = double(specific) ; convert specific to watervapor pressure vapor_pressure=dblarr(size(specific, /dim)) ; generate water vapor pressure, for refractivity for pressure_loop = 0, n_elements(pressure)-1 do begin vapor_pressure[*,*,pressure_loop] = pressure[pressure_loop]*specific[*,*,pressure_loop]/(0.622+0.378*specific[*,*,pressure_loop]) endfor ; convert kg/kg to g/kg specific *= 1000.D END FUNCTION SUB_INPUT_CFS, pdatestr_dashed, $ ; ------- lon=lon_cfs, $ lat=lat_cfs, $ pressure=pressure_cfs, $ ; ------- Geopotential_height_all, $ Temperature_all, $ watervapor_all, $ ref=refractivity_all, $ dryt=dryTemperature_all, $ type = type ; h06/h00 ;print,'pdatestr_dashed',pdatestr_dashed pdatestr_array=strsplit(pdatestr_dashed,'-',/extract) pyearstr=pdatestr_array[0] pmonthstr=pdatestr_array[1] pdaystr=pdatestr_array[2] pdatestr = pyearstr+pmonthstr+pdaystr caldat,julday(pmonthstr,fix(pdaystr)+1,pyearstr), pmonthstr_tomorrow, pdaystr_tomorrow, pyearstr_tomorrow pyearstr_tomorrow = string(pyearstr_tomorrow, format='(i4)') pdatestr_tomorrow = pyearstr_tomorrow+string(pmonthstr_tomorrow, format='(i2.2)')+string(pdaystr_tomorrow, format='(i2.2)') pdatestr_tomorrow_dashed = pyearstr_tomorrow+'-'+string(pmonthstr_tomorrow, format='(i2.2)')+'-'+string(pdaystr_tomorrow, format='(i2.2)') caldat,julday(pmonthstr,fix(pdaystr)-1,pyearstr), pmonthstr_yesterday, pdaystr_yesterday, pyearstr_yesterday pyearstr_yesterday = string(pyearstr_yesterday, format='(i4)') pdatestr_yesterday = pyearstr_yesterday+string(pmonthstr_yesterday, format='(i2.2)')+string(pdaystr_yesterday, format='(i2.2)') pdatestr_yesterday_dashed = pyearstr_yesterday+'-'+string(pmonthstr_yesterday, format='(i2.2)')+'-'+string(pdaystr_yesterday, format='(i2.2)') ; ---------------------------- spawn, 'hostname -s', hostname if strmid(hostname, 0, 3) eq 'rhw' then begin folder= '/data/data504/xzhou/' endif else if hostname eq 'gnssro' then begin folder= '/data/xinjiaz/' endif else if hostname eq 'starro' then begin folder= '/GNSSRO_data/xinjiaz/' endif else if hostname eq 'gpsmet' then begin folder= '/GNSSRO_data/xinjiaz/' endif L4_patten = strarr(9) if pdatestr_yesterday_dashed le '2011-03-31' then begin L4_patten[0] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/pgb'+type+'.gdas.'+pdatestr_yesterday+'00.nc' L4_patten[1] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/pgb'+type+'.gdas.'+pdatestr_yesterday+'06.nc' L4_patten[2] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/pgb'+type+'.gdas.'+pdatestr_yesterday+'12.nc' L4_patten[3] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/pgb'+type+'.gdas.'+pdatestr_yesterday+'18.nc' endif else begin L4_patten[0] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/cdas1.t00z.pgrb'+type+'.nc' L4_patten[1] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/cdas1.t06z.pgrb'+type+'.nc' L4_patten[2] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/cdas1.t12z.pgrb'+type+'.nc' L4_patten[3] = folder+'data_CFS/'+type+'/'+pdatestr_yesterday_dashed+'/cdas1.t18z.pgrb'+type+'.nc' endelse if pdatestr_dashed le '2011-03-31' then begin L4_patten[4] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/pgb'+type+'.gdas.'+pdatestr+'00.nc' L4_patten[5] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/pgb'+type+'.gdas.'+pdatestr+'06.nc' L4_patten[6] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/pgb'+type+'.gdas.'+pdatestr+'12.nc' L4_patten[7] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/pgb'+type+'.gdas.'+pdatestr+'18.nc' endif else begin L4_patten[4] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/cdas1.t00z.pgrb'+type+'.nc' L4_patten[5] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/cdas1.t06z.pgrb'+type+'.nc' L4_patten[6] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/cdas1.t12z.pgrb'+type+'.nc' L4_patten[7] = folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/cdas1.t18z.pgrb'+type+'.nc' endelse if pdatestr_tomorrow_dashed le '2011-03-31' then begin L4_patten[8] = folder+'data_CFS/'+type+'/'+pdatestr_tomorrow_dashed+'/pgb'+type+'.gdas.'+pdatestr_tomorrow+'00.nc' endif else begin L4_patten[8] = folder+'data_CFS/'+type+'/'+pdatestr_tomorrow_dashed+'/cdas1.t00z.pgrb'+type+'.nc' endelse if type eq 'h24' then L4_patten_this = L4_patten[0:4] $ else if type eq 'h18' then L4_patten_this = L4_patten[1:5] $ else if type eq 'h12' then L4_patten_this = L4_patten[2:6] $ else if type eq 'h06' then L4_patten_this = L4_patten[3:7] $ ; < --- use 6 hour forecast else if type eq 'h00' then L4_patten_this = L4_patten[4:8] $ else stop ; for 'f003',f009',f015',f021', the time_loc is not integer, so do not use these 4 forecast_interval ; L4_FileList = FILE_SEARCH(L4_patten_this, COUNT=L4_patten_Count) print, '' if L4_patten_Count ne 5 then begin print, 'there is no enough CFS('+type+') files in '+pdatestr_dashed print, folder+'data_CFS/'+type+'/'+pdatestr_dashed+'/cdas1.t??z.pgrb'+type+'.nc' print, folder+'data_CFS/'+type+'/'+pdatestr_tomorrow_dashed+'/cdas1.t00z.pgrb'+type+'.nc' ;print, L4_patten_this ;print, '----------' temperature_all=!VALUES.F_NAN watervapor_all=!VALUES.F_NAN Geopotential_height_all=!VALUES.F_NAN dryTemperature_all=!VALUES.F_NAN refractivity_all=!VALUES.F_NAN return, 0 endif else begin restore,folder+'data_CFS/dim.sav' ;lon_cfs, [ 0 ~ 360], 720 ;lat_cfs, [90 ~ -90], 361 ;pressure_cfs, [1 ~ 1000], 37 levels, hPa temperature_all= [] watervapor_all = [] vapor_pressure_all = [] Geopotential_height_all=[] for file_loop = 0, L4_patten_Count-1 do begin ; -------------------- print,L4_FileList[file_loop] Read_reference_CFS, L4_FileList[file_loop], $ pressure_cfs, $ temperature, $ watervapor, $ vapor_pressure, $ Geopotential_height temperature = [temperature, temperature[0, *, *]] ; add one more column for longitude watervapor = [watervapor, watervapor[0, *, *]] vapor_pressure = [vapor_pressure, vapor_pressure[0, *, *]] Geopotential_height = [Geopotential_height, Geopotential_height[0, *, *]] size_dim = size(temperature, /dim) temperature_all = [temperature_all, reform(temperature, [1, size_dim[0], size_dim[1], size_dim[2]])] watervapor_all = [watervapor_all, reform(watervapor, [1, size_dim[0], size_dim[1], size_dim[2]])] vapor_pressure_all = [vapor_pressure_all, reform(vapor_pressure, [1, size_dim[0], size_dim[1], size_dim[2]])] Geopotential_height_all= [Geopotential_height_all, reform(Geopotential_height, [1, size_dim[0], size_dim[1], size_dim[2]])] endfor ; for interpolation, mv 'vertical' to left most, because interpolation starting from rightmost Temperature_all = transpose(temperature_all, [3, 0, 1, 2]) ; 31:5:360:181, height:time:lon:lat watervapor_all = transpose(watervapor_all, [3, 0, 1, 2]) vapor_pressure_all = transpose(vapor_pressure_all, [3, 0, 1, 2]) Geopotential_height_all = transpose(Geopotential_height_all, [3, 0, 1, 2]) ; rebuild pressure size_dim = size(Temperature_all, /dim) pressure_all = dblarr(size_dim) for ii = 0, size_dim[0]-1 do pressure_all[ii,*,*,*]=pressure_cfs[ii] ; calculate refractivity if ARG_PRESENT(refractivity_all) then refractivity_all = 77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all ; calculate dry temperature if ARG_PRESENT(dryTemperature_all) then dryTemperature_all = 77.6D*pressure_all/(77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all) return, 2 endelse END