PRO Read_reference_JRA, input_file, $ pdaystr_this, $ pressure, $ temperature, $ watervapor, $ ; specific humidity vapor_pressure, $ ; for refractivity Geopotential_height ncfileid = NCDF_OPEN(input_file) print, 'TMP_GDS0_ISBL:',input_file NCDF_VARGET, ncfileid, 'TMP_GDS0_ISBL', temperature ; K, [288, 145, 37, 124], float NCDF_CLOSE, ncfileid input_file = StrJoin(StrSplit(input_file,'011_tmp',/Regex,/Extract,/Preserve_Null), '051_spfh') ncfileid = NCDF_OPEN(input_file) print, 'SPFH_GDS0_ISBL:',input_file NCDF_VARGET, ncfileid, 'SPFH_GDS0_ISBL', specific_27; kg/kg, [288, 145, 27, 124], float NCDF_CLOSE, ncfileid size_t = size(temperature, /dim) if pdaystr_this eq 0 then begin this_id = 0 ; read L4_patten_tomorrow, pick only the first slot endif else begin if pdaystr_this*4 eq size_t[3] then this_id = [(pdaystr_this-1)*4:pdaystr_this*4-1 ] $ ; the last day, only pick 4 times else this_id = [(pdaystr_this-1)*4:pdaystr_this*4] ; not the last day, pick 5 times endelse ;print,'size(temperature), inside 1:',size(temperature) temperature = double(temperature[*,*, *, this_id]) ;print,'size(temperature), inside 2:',size(temperature) specific = dblarr(288, 145, 37, n_elements(this_id)) valid_id_37 = indgen(27)+10 ; specific_27 is 10 levels less than temperature, which is the top 10 levels specific[*, *, valid_id_37, *] = specific_27[*, *, *, this_id] watervapor=specific*1000. vapor_pressure=[specific] 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 input_file = StrJoin(StrSplit(input_file,'051_spfh',/Regex,/Extract,/Preserve_Null), '007_hgt') ncfileid = NCDF_OPEN(input_file) print, 'HGT_GDS0_ISBL:',input_file NCDF_VARGET, ncfileid, 'HGT_GDS0_ISBL', Geopotential_height ; gpm, [288, 145, 37, 124], float NCDF_CLOSE, ncfileid Geopotential_height = Geopotential_height[*,*, *, this_id] END FUNCTION SUB_INPUT_JRA, pdatestr_dashed, $ ; ------- lon=lon_jra, $ lat=lat_jra, $ pressure=pressure_jra, $ ; ------- Geopotential_height_all, $ Temperature_all, $ watervapor_all, $ ref=refractivity_all, $ dryt=dryTemperature_all ;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 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)') ; ---------------------------- 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 ; 011_tmp, 007_ght, 051_spfh L4_patten_today = folder+'data_JRA/data/'+pyearstr+ '/anl_p125.011_tmp.'+pdatestr+ '0100_'+pdatestr+ '??18.nc' L4_patten_tomorrow = folder+'data_JRA/data/'+pyearstr_tomorrow+'/anl_p125.011_tmp.'+pdatestr_tomorrow+'0100_'+pdatestr_tomorrow+'??18.nc' L4_patten = [L4_patten_today, L4_patten_tomorrow] L4_patten = L4_patten[uniq(L4_patten)] L4_FileList = FILE_SEARCH(L4_patten, COUNT=L4_patten_Count) print, '' if L4_patten_Count eq 0 then begin print, 'there is no enough JRA files in '+pdatestr_dashed+' '+L4_patten_today 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_JRA/dim.sav' ;lon_jra, [ 0 ~ 360] ;lat_jra, [90 ~ -90] ;pressure_jra, [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] if file_loop eq 0 then pdaystr_this = fix(pdaystr) $ ; for reading L4_patten_today else pdaystr_this = 0 ; for reading L4_patten_tomorrow ;print, 'pdaystr_this=',pdaystr_this Read_reference_JRA, L4_FileList[file_loop] ,pdaystr_this, pressure_jra, 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, *, *, *]] ;print,'size(temperature)', size(temperature) ;print,'size(watervapor)', size(watervapor) if file_loop eq 0 then begin ; normal mode, 4 dimension temperature_all = [temperature_all, transpose(temperature, [3, 0, 1, 2])] ; put time in first dim, for combining watervapor_all = [watervapor_all, transpose(watervapor, [3, 0, 1, 2])] ; now it's [4, 288, 145, 37] vapor_pressure_all = [vapor_pressure_all, transpose(vapor_pressure, [3, 0, 1, 2])] ; now it's [4, 288, 145, 37] Geopotential_height_all= [Geopotential_height_all, transpose(Geopotential_height, [3, 0, 1, 2])] endif else begin ; next month, 3 dimension temperature_all = [temperature_all, reform(temperature, [1, 289, 145, 37])] ; put time in first dim, for combining watervapor_all = [watervapor_all, reform(watervapor, [1, 289, 145, 37])] ; now it's [1, 288, 145, 37] vapor_pressure_all = [vapor_pressure_all, reform(vapor_pressure, [1, 289, 145, 37])] ; now it's [1, 288, 145, 37] Geopotential_height_all= [Geopotential_height_all, reform(Geopotential_height, [1, 289, 145, 37])] endelse endfor ; for interpolation, mv 'vertical' to left most, because interpolation starting from rightmost Temperature_all = transpose(temperature_all, [3, 0, 1, 2]) ; 37:5:288:145, 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 for refractivity and dry temperature size_dim = size(Temperature_all, /dim) pressure_all = dblarr(size_dim) for ii = 0, size_dim[0]-1 do pressure_all[ii,*,*,*]=pressure_jra[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