PRO Read_reference_GFS, input_file, $ pressure, $ temperature, $ watervapor, $ ; specific humidity vapor_pressure, $ ; for refractivity Geopotential_height ncfileid = NCDF_OPEN(input_file) temperature_id1=NCDF_VARID(ncfileid, 'TMP_P0_L100_GLL0') if temperature_id1 gt 0 then begin NCDF_VARGET, ncfileid, 'TMP_P0_L100_GLL0', temperature ; K, [360, 181, 31], float NCDF_VARGET, ncfileid, 'RH_P0_L100_GLL0', RelativeHumidity ; %, [360, 181, 31], float NCDF_VARGET, ncfileid, 'HGT_P0_L100_GLL0', Geopotential_height ; gpm, [360, 181, 17], float endif else begin NCDF_VARGET, ncfileid, 'TMP_3_ISBL', temperature ; K, [360, 181, 31], float NCDF_VARGET, ncfileid, 'R_H_3_ISBL', RelativeHumidity ; %, [360, 181, 31], float NCDF_VARGET, ncfileid, 'HGT_3_ISBL', Geopotential_height ; gpm, [360, 181, 17], float endelse NCDF_CLOSE, ncfileid ; GFS format, temperature/Geopotential_height/RelativeHumidity ; 2021-03-22 ~ current : (lv_ISBL0 = 41) ; 2019-06-13 ~ 2021-03-22: (lv_ISBL0 = 34 / lv_ISBL5 = 31) ** 2021-03-22 06:00|12:00 ; 2016-05-12 ~ 2019-06-12: (lv_ISBL0 = 31) <---------------- default ; 2015-03-24 ~ 2016-05-11: (lv_ISBL0 = 26) ; 2015-01-15 ~ 2015-03-23: missing ; 2012-05-22 ~ 2015-01-14: (lv_ISBL0 = 26 / lv_ISBL6 = 25) ** 2012-05-22 06:00|12:00 ; 2009-12-15 ~ 2012-05-22: (lv_ISBL0 = 26 / lv_ISBL4 = 25) ; 2008-02-01 ~ 2009-12-15: (lv_ISBL0 = 26 / lv_ISBL4 = 21) ** 2009-12-15 06:00|12:00 ; 2004-03-02 ~ 2008-01-31: (lv_ISBL3 = 26 / lv_ISBL7 = 21) size_t=size(temperature, /dim) size_wv=size(RelativeHumidity, /dim) if size_t[2] eq 41 then begin ; (lv_ISBL0 = 41) ; in this format, temperature/geopotential/RelativeHumidity extended to 41 levels valid_id = [[8:13],15,16,[18:40]] temperature = temperature[*,*,valid_id] Geopotential_height = Geopotential_height[*,*,valid_id] RelativeHumidity = RelativeHumidity[*,*,valid_id] endif else if size_t[2] eq 34 then begin ; (lv_ISBL0 = 34 / lv_ISBL5 = 31) ; in this format, temperature/geopotential extended to 34 levels, RelativeHumidity remain 31 levels valid_id=[1,2,3,4,5,6,8,9,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33] temperature = temperature[*,*,valid_id] Geopotential_height = Geopotential_height[*,*,valid_id] endif else if size_t[2] eq 26 and size_wv[2] eq 26 then begin ; (lv_ISBL0 = 26) ; in this format, temperature/RelativeHumidity/geopotential are in 26 levels, 10hPh to 1000hPa valid_id_26=indgen(26)+5 temperature_26 = temporary(temperature) RelativeHumidity_26= temporary(RelativeHumidity) Geopotential_height_26= temporary(Geopotential_height) temperature = fltarr(360, 181, 31)+!values.f_nan RelativeHumidity = fltarr(360, 181, 31)+!values.f_nan Geopotential_height= fltarr(360, 181, 31)+!values.f_nan ;print, input_file ;print, 'size(temperature_26)' ;print, size(temperature_26) temperature[*,*,valid_id_26] = temperature_26 RelativeHumidity[*,*,valid_id_26] = RelativeHumidity_26 Geopotential_height[*,*,valid_id_26] = Geopotential_height_26 endif else if size_t[2] eq 26 and size_wv[2] eq 25 then begin ; (lv_ISBL0 = 26 / lv_ISBL6 = 25), or (lv_ISBL0 = 26 / lv_ISBL4 = 25) ; in this format, temperature/RelativeHumidity/geopotential are in 26 levels, 10hPh to 1000hPa valid_id_26=indgen(26)+5 valid_id_25=indgen(25)+6 & valid_id_25[0]=5 temperature_26 = temporary(temperature) RelativeHumidity_25= temporary(RelativeHumidity) Geopotential_height_26= temporary(Geopotential_height) temperature = fltarr(360, 181, 31)+!values.f_nan RelativeHumidity = fltarr(360, 181, 31) RelativeHumidity[*, *, valid_id_25] = !values.f_nan ; keep water in 20hPa as zero Geopotential_height= fltarr(360, 181, 31)+!values.f_nan temperature[*,*,valid_id_26] = temperature_26 RelativeHumidity[*,*,valid_id_25] = RelativeHumidity_25 Geopotential_height[*,*,valid_id_26] = Geopotential_height_26 endif else if size_t[2] eq 26 and size_wv[2] eq 21 then begin ; (lv_ISBL0 = 26 / lv_ISBL4 = 21), or (lv_ISBL3 = 26 / lv_ISBL7 = 21) ; in this format, temperature/RelativeHumidity/geopotential are in 26 levels, 10hPh to 1000hPa valid_id_26=indgen(26)+5 valid_id_21=indgen(21)+10 temperature_26 = temporary(temperature) RelativeHumidity_21= temporary(RelativeHumidity) Geopotential_height_26= temporary(Geopotential_height) temperature = fltarr(360, 181, 31)+!values.f_nan RelativeHumidity = fltarr(360, 181, 31) RelativeHumidity[*, *, valid_id_21] = !values.f_nan ; keep water in 20hPa as zero Geopotential_height= fltarr(360, 181, 31)+!values.f_nan ;print, input_file ;print, 'size(temperature_26)' ;print, size(temperature_26) temperature[*,*,valid_id_26] = temperature_26 RelativeHumidity[*,*,valid_id_21] = RelativeHumidity_21 Geopotential_height[*,*,valid_id_26] = Geopotential_height_26 endif size_r=size(RelativeHumidity, /dim) if size_r[2] eq 25 then begin ; GDAS from theia format, 25 levels valid_id_31=[5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30] RelativeHumidity_25 = RelativeHumidity RelativeHumidity = fltarr(360, 181, 31)+!values.f_nan RelativeHumidity[*,*,valid_id_31]= RelativeHumidity_25 endif temperature = double(temperature) saturated_vp = 6.11 * (10^ (7.5 * (temperature-273.15) / (237.7 + (temperature-273.15)))) vapor_pressure = RelativeHumidity / 100.D * saturated_vp WaterVapor = fltarr(360, 181, 31)+!values.f_nan for pressure_loop = 0, n_elements(pressure)-1 do begin WaterVapor[*,*,pressure_loop] = 622.*Vapor_pressure[*,*,pressure_loop]/(pressure[pressure_loop] - 0.378*Vapor_pressure[*,*,pressure_loop]) ; g/kg endfor END FUNCTION SUB_INPUT_GFS, pdatestr_dashed, $ ; ------- lon=lon_gfs, $ lat=lat_gfs, $ pressure=pressure_gfs, $ ; ------- Geopotential_height_all, $ Temperature_all, $ watervapor_all, $ ref=refractivity_all, $ dryt=dryTemperature_all, $ type = type ; GFS/GDAS ;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)') 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)') if keyword_set(type) then begin if type eq 'GFS' then begin folder_name = 'GFS_h06' type2 = 'gfs' forecast_interval = 'f006' ; ['f000','f003','f006','f009','f012','f015','f018','f021','f024'] endif else if type eq 'GDAS' then begin folder_name = 'GFS_h00' type2 = '*' forecast_interval = 'f000' endif endif else begin ; default is GFS type = 'GFS' folder_name = 'GFS_h06' type2 = 'gfs' forecast_interval = 'f006' ; ['f000','f003','f006','f009','f012','f015','f018','f021','f024'] endelse ; ---------------------------- 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) L4_patten[0] = folder+'data_GFS/'+folder_name+'/'+pyearstr_yesterday+'/'+type2+'.t00z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr_yesterday+'.nc' L4_patten[1] = folder+'data_GFS/'+folder_name+'/'+pyearstr_yesterday+'/'+type2+'.t06z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr_yesterday+'.nc' L4_patten[2] = folder+'data_GFS/'+folder_name+'/'+pyearstr_yesterday+'/'+type2+'.t12z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr_yesterday+'.nc' L4_patten[3] = folder+'data_GFS/'+folder_name+'/'+pyearstr_yesterday+'/'+type2+'.t18z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr_yesterday+'.nc' L4_patten[4] = folder+'data_GFS/'+folder_name+'/'+pyearstr+'/'+type2+'.t00z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr+'.nc' L4_patten[5] = folder+'data_GFS/'+folder_name+'/'+pyearstr+'/'+type2+'.t06z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr+'.nc' L4_patten[6] = folder+'data_GFS/'+folder_name+'/'+pyearstr+'/'+type2+'.t12z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr+'.nc' L4_patten[7] = folder+'data_GFS/'+folder_name+'/'+pyearstr+'/'+type2+'.t18z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr+'.nc' L4_patten[8] = folder+'data_GFS/'+folder_name+'/'+pyearstr_tomorrow+'/'+type2+'.t00z.pgrb2.1p00.'+forecast_interval+'.'+pdatestr_tomorrow+'.nc' ; generating_time forecast_interval ; forecast_time = generating_time + forecast_interval if forecast_interval eq 'f024' then L4_patten_this = L4_patten[0:4] $ else if forecast_interval eq 'f018' then L4_patten_this = L4_patten[1:5] $ else if forecast_interval eq 'f012' then L4_patten_this = L4_patten[2:6] $ else if forecast_interval eq 'f006' then L4_patten_this = L4_patten[3:7] $ ; < --- use 6 hour forecast else if forecast_interval eq 'f000' 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 GFS('+type+') files in '+pdatestr_dashed 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_GFS/dim_1P00.sav' ;lon_gfs, [ 0 ~ 360], 360 ;lat_gfs, [90 ~ -90], 181 ;pressure_gfs, [1 ~ 1000], 31 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_GFS, L4_FileList[file_loop], pressure_gfs, 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, *, *]] temperature_all = [temperature_all, reform(temperature, [1, 361, 181, 31])] watervapor_all = [watervapor_all, reform(watervapor, [1, 361, 181, 31])] vapor_pressure_all = [vapor_pressure_all, reform(vapor_pressure, [1, 361, 181, 31])] Geopotential_height_all= [Geopotential_height_all, reform(Geopotential_height, [1, 361, 181, 31])] 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_gfs[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