PRO Read_reference_MERRA, input_file, $ pressure, $ temperature, $ watervapor, $ ; specific humidity vapor_pressure, $ ; for refractivity Geopotential_height, $ for_tomorrow=for_tomorrow ncfileid = NCDF_OPEN(input_file) NCDF_VARGET, ncfileid, 'T', temperature ; K, [576, 361, 42, 8], double NCDF_VARGET, ncfileid, 'QV', specific ; kg/kg, [576, 361, 42, 8], double NCDF_VARGET, ncfileid, 'H', Geopotential_height ; gpm, [144, 73, 17], float ;NCDF_VARGET, ncfileid, 'lev', pressure_merra ; 42, hPa [1000, 0.1] ;pressure_merra = reverse(float(pressure_merra)) ; revers, -> [0.1, 1000] ;NCDF_VARGET, ncfileid, 'lat', lat_merra ; 361, [ -90 ~ 90] ;NCDF_VARGET, ncfileid, 'lon', lon_merra ; 576, [-180 ~ 180] NCDF_CLOSE, ncfileid if keyword_set(for_tomorrow) then begin ; only keep the first time temperature = temperature[*, *, *, 0] specific = specific[*, *, *, 0] Geopotential_height = Geopotential_height[*, *, *, 0] endif temperature = reverse(temperature, 3) ; reverse the pressure, specific = reverse(specific, 3) ; [1000 ~ 0.1] -> [0.1 ~ 1000] Geopotential_height = reverse(Geopotential_height, 3) ; [1000 ~ 0.1] -> [0.1 ~ 1000] ; ----------------- valid_id=where(temperature gt 9.e14, valid_num) if valid_num gt 0 then temperature[valid_id] = !values.d_nan valid_id=where(specific gt 9.e14, valid_num) if valid_num gt 0 then specific[valid_id] = !values.d_nan valid_id=where(Geopotential_height gt 9.e14, valid_num) if valid_num gt 0 then Geopotential_height[valid_id] = !values.d_nan ; ----------------- ; convert kg/kg to g/kg watervapor=specific*1000. ; convert specific to watervapor pressure 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 END FUNCTION SUB_INPUT_MERRA, pdatestr_dashed, $ ; ------- lon=lon_merra, $ lat=lat_merra, $ pressure=pressure_merra, $ ; ------- 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+pdaystr caldat,julday(pmonthstr,fix(pdaystr)+1,pyearstr), pmonthstr_tomorrow, pdaystr_tomorrow, pyearstr_tomorrow pyearstr_tomorrow = string(pyearstr_tomorrow, format='(i4)') pmonthstr_tomorrow = string(pmonthstr_tomorrow, format='(i2.2)') pdatestr_tomorrow = pyearstr_tomorrow+pmonthstr_tomorrow+string(pdaystr_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= '/STARRO_data2/xinjiaz/' endif else if hostname eq 'starro' then begin folder= '/data2/xinjiaz/' endif else if hostname eq 'gpsmet' then begin folder= '/STARRO_data2/xinjiaz/' endif L4_file_today = file_search(folder+'data_MERRA/data/M2I3NPASM.5.12.4/'+pyearstr +'/'+pmonthstr +'/MERRA2_?0?.inst3_3d_asm_Np.'+pdatestr +'.nc4', count=count_today) L4_file_tomorrow = file_search(folder+'data_MERRA/data/M2I3NPASM.5.12.4/'+pyearstr_tomorrow+'/'+pmonthstr_tomorrow+'/MERRA2_?0?.inst3_3d_asm_Np.'+pdatestr_tomorrow+'.nc4', count=count_tomorrow) print, '' if count_today eq 0 or count_tomorrow eq 0 then begin print, 'there is no enough MERRA files in: '+folder+'data_MERRA/data/M2I3NPASM.5.12.4/'+pyearstr +'/'+pmonthstr +'/MERRA2_?0?.inst3_3d_asm_Np.'+pdatestr +'.nc4' print, 'there is no enough MERRA files in: '+folder+'data_MERRA/data/M2I3NPASM.5.12.4/'+pyearstr_tomorrow+'/'+pmonthstr_tomorrow+'/MERRA2_?0?.inst3_3d_asm_Np.'+pdatestr_tomorrow+'.nc4' 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_MERRA/dim.sav' ;lon_merra 576+1, [-180 ~ 180] ;lat_merra, 361, [ -90 ~ 90] ;pressure_merra, 42, [0.1 ~ 1000 ] temperature_all= [] watervapor_all = [] vapor_pressure_all = [] Geopotential_height_all=[] ; -------------------- print, L4_file_today Read_reference_MERRA, L4_file_today[0], $ pressure_merra, $ temperature, $ watervapor, $ vapor_pressure, $ Geopotential_height temperature = [temperature, temperature[0, *, *, *]] ; add one more column for longitude watervapor = [watervapor, watervapor[0, *, *, *]] ; [576, 361, 42, 8] -> [577, 361, 42, 8] vapor_pressure = [vapor_pressure, vapor_pressure[0, *, *, *]] ; [576, 361, 42, 8] -> [577, 361, 42, 8] Geopotential_height = [Geopotential_height, Geopotential_height[0, *, *, *]] Temperature = transpose(temperature, [3, 0, 1, 2]) ; mv time to first column, for combining with data in tomorrow watervapor = transpose(watervapor, [3, 0, 1, 2]) ; [577, 361, 42, 8] -> [8, 577, 361, 42] vapor_pressure = transpose(vapor_pressure, [3, 0, 1, 2]) ; [577, 361, 42, 8] -> [8, 577, 361, 42] Geopotential_height = transpose(Geopotential_height, [3, 0, 1, 2]) ; -------------------- print, L4_file_tomorrow Read_reference_MERRA, L4_file_tomorrow[0], $ pressure_merra, $ temperature_tomorrow, $ watervapor_tomorrow, $ vapor_pressure_tomorrow, $ Geopotential_height_tomorrow, $ for_tomorrow=1 temperature_tomorrow = [temperature_tomorrow, temperature_tomorrow[0, *, *]] ; add one more column for longitude watervapor_tomorrow = [watervapor_tomorrow, watervapor_tomorrow[0, *, *]] ; [576, 361, 42] -> [577, 361, 42] vapor_pressure_tomorrow = [vapor_pressure_tomorrow, vapor_pressure_tomorrow[0, *, *]] ; [576, 361, 42] -> [577, 361, 42] Geopotential_height_tomorrow = [Geopotential_height_tomorrow, Geopotential_height_tomorrow[0, *, *]] temperature_tomorrow = reform(temperature_tomorrow, [1, 577, 361, 42]) ; add one column, for combining data watervapor_tomorrow = reform(watervapor_tomorrow, [1, 577, 361, 42]) vapor_pressure_tomorrow = reform(vapor_pressure_tomorrow, [1, 577, 361, 42]) Geopotential_height_tomorrow = reform(Geopotential_height_tomorrow, [1, 577, 361, 42]) ; -------------------- temperature_all= [Temperature, temperature_tomorrow] ; combine today and tomorrow watervapor_all = [watervapor, watervapor_tomorrow] ; [8, 577, 361, 42] -> [9, 577, 361, 42] vapor_pressure_all = [vapor_pressure, vapor_pressure_tomorrow] ; [8, 577, 361, 42] -> [9, 577, 361, 42] Geopotential_height_all= [Geopotential_height, Geopotential_height_tomorrow] ; combine today and tomorrow ; -------------------- Temperature_all = transpose(temperature_all, [3, 0, 1, 2]) ; mv 'vertical' to left most, because interpolation starting from rightmost watervapor_all = transpose(watervapor_all, [3, 0, 1, 2]) ; [9, 577, 361, 42] -> [42, 9, 577, 361] vapor_pressure_all = transpose(vapor_pressure_all, [3, 0, 1, 2]) ; [9, 577, 361, 42] -> [42, 9, 577, 361] Geopotential_height_all = transpose(Geopotential_height_all, [3, 0, 1, 2]) ; mv 'vertical' to left most, because interpolation starting from rightmost ; ------------------- ; 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_merra[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.6*pressure_all/(77.6D*pressure_all/Temperature_all + 3.73e5 * vapor_pressure_all / Temperature_all / Temperature_all) return, 2 endelse END