FUNCTION Extract_UCARbufr, input_file , $ data_this_file, $ version=version return_code=0B ; 2 is good ; 1 is bad but no email ; 0 is bad and send email if keyword_set(version) then begin mission=version endif else begin path=strsplit(file_dirname(input_file),'/',/extract) mission=path[-4] endelse agency = 'UCARbufr' filename_array=strsplit(file_basename(input_file), '_.', /extract) leo_name = filename_array[1] & if mission eq 'spire' then leo_name = strmid(leo_name, 1) ; "S102" --> "102" gps_name = strlowcase(filename_array[6]) data_this_file=[] ; make sure the file name format is correct ('bfrPrf_C006.2018.305.03.54.G10_2016.1120') if strmid(file_basename(input_file), 26, 1) ne '.' or $ strmid(file_basename(input_file), 30, 1) ne '_' then begin error_string=input_file+' is in different filename format' goto, jump_end endif ;print,input_file cmd='openr,11,input_file' res=execute(cmd) if res eq 0 then begin error_string=input_file+' cannot be openned' goto, jump_end endif ; COMMON CODE TABLE C–5 ; equavelant to 'cic085' leo_list = [['MetopB', '3.0' ],$ ['MetopA', '4.0' ],$ ['MetopC', '5.0' ],$ ['CHAMP', '41.0' ],$ ['TerraSAR-X', '42.0' ],$ ['TanDem-X', '43.0' ],$ ['PAZ', '44.0' ],$ ['GRACE-A', '722.0' ],$ ['GRACE-B', '723.0' ],$ ['COSMIC-1', '740.0'],$ ['COSMIC-2', '741.0'],$ ['COSMIC-3', '742.0'],$ ['COSMIC-4', '743.0'],$ ['COSMIC-5', '744.0'],$ ['COSMIC-6', '745.0'],$ ['C2E1', '750.0'],$ ['C2E2', '751.0'],$ ['C2E3', '752.0'],$ ['C2E4', '753.0'],$ ['C2E5', '754.0'],$ ['C2E6', '755.0'],$ ['KOMPSAT5', '825.0'],$ ['C-NOFS', '786.0' ]] ; equal to 'g01'/'r01' gps_list = [['GPS', '401.0' ],$ ['GLONASS', '402.0' ],$ ['GALILEO', '403.0' ],$ ['BDS', '404.0' ]] line_num = file_lines(input_file) file_content = strarr(line_num) readf,11,file_content close,11 end_id=where(strpos(file_content, 'Reached end of BUFR file') ge 0) ; only one file_content = file_content[0 : end_id[0]-1] ; --- title (cic085_2018-11-30T08:16:00_g01_GeoOptics_CICERO) leo_id=where(strpos(file_content, '001007') eq 0) ;MetopA ;agency_instrument_id=where(strpos(file_content, '002019') eq 0) ; ESA/EUMETSAT; GRAS; GNSS receiver gps_id=where(strpos(file_content, '002020') eq 0) ; GPS/GLONAS/GALILIO/BDS QFRO_id=where(strpos(file_content, '033039') eq 0) ; see https://www.nco.ncep.noaa.gov/sib/jeff/CodeFlag_0_STDv14_LOC7.html # 033039 , 0-33-039 - QFRO ELRC_id=where(strpos(file_content, '010035') eq 0) ; Earth's local radius of curvature, unit is meter year_id=where(strpos(file_content, '004001') eq 0) month_id=where(strpos(file_content, '004002') eq 0) day_id=where(strpos(file_content, '004003') eq 0) hour_id=where(strpos(file_content, '004004') eq 0) minute_id=where(strpos(file_content, '004005') eq 0) second_id=where(strpos(file_content, '004006') eq 0) ; --- leo leo_str = strsplit(file_content[leo_id[0]], /extract) ; 001007 ;print,'file_content[leo_id[0]]: ', file_content[leo_id[0]] ;print,'leo_str: ',leo_str if leo_str[2] eq '999.0' or mission eq 'cicero' or mission eq 'spire' or mission eq 'planetiqrt' then begin leo = leo_name endif else begin this_id = where(leo_list[1,*] eq leo_str[2], this_num) if this_num eq 0 then begin cmd='echo -e "Search '+leo_str[2]+' ('+file_content[leo_id[0]]+') at "' ;cmd=cmd+'" \n COMMON CODE TABLE C–5 at https://www.wmo.int/pages/prog/www/WMOCodes/WMO306_vI2/LatestVERSION/WMO306_vI2_CommonTable_en.pdf"' cmd=cmd+'" \n COMMON CODE TABLE C–5 at https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/WMO306_vI2_CommonTable_en_v23.0.0.pdf"' cmd=cmd+'" \n decode file: '+input_file+'"' cmd=cmd+' | mail -s "Alert: New leo satellite in Extract_UCARbufr.pro " xinjia.zhou@noaa.gov' spawn, cmd stop endif leo = leo_list[0,this_id[0]] endelse ; --- gps if !false then begin ; --- read gps name from bufr table gps_str = strsplit(file_content[gps_id[0]], /extract) this_id = where(gps_list[1,*] eq gps_str[2], this_num) if this_num eq 0 then begin cmd=string('echo -e "Search --'+gps_str+'-- at \n Satellite classification at https://www.wmo.int/pages/prog/www/WMOCodes/WMO306_vI2/LatestVERSION/WMO306_vI2_BUFRCREX_CodeFlag_en.pdf" | mail -s "Alert: New gps in Extract_UCARbufr.pro " xinjia.zhou@noaa.gov') spawn, cmd endif gps = gps_list[0,this_id[0]] endif else begin ; --- read gps name from filename gps = gps_name endelse ; --- bad QFRO_str = file_content[QFRO_id[0]] ; 033039 QFRO_array = strsplit(QFRO_str, /extract) QFRO_value = QFRO_array[2] ; example: 256.0 -> [8] ; 8448.0 -> [3,8] ; 8512.0 -> [3,8,10] ; 61184.0 -> [1,2,3,5,6,7,8] QFRO_bit=[] for ii = 15, 0,-1 do begin ; the sequence is reversed in WMO bufr if 2.^ii gt QFRO_value then begin continue endif else begin QFRO_bit = [QFRO_bit, 16-ii] QFRO_value -= 2.^ii endelse endfor if total(QFRO_bit eq 1 or $ ; Non-nominal quality QFRO_bit eq 4 or $ ; Excess Phase processing non-nominal QFRO_bit eq 5 or $ ; Bending Angle processing non-nominal QFRO_bit eq 6 or $ ; Refractivity processing non-nominal QFRO_bit eq 7 or $ ; Meteorological processing non-nominal QFRO_bit eq 14) $ ; Background profile non-nominal gt 0 then is_bad = 1 else is_bad = 0 ;is_bad = strmatch(QFRO_str, '*MANY BITS ON*') ; this one is not always working, see example at: cosmic2013/bfrPrf/decode/2012.001/bfrPrf_C001.2012.001.00.06.G28_2013.3520 if total(QFRO_bit eq 3) gt 0 then irs = -1 else irs = 1 ; Ascending occultation flag ; --- time year_str = strsplit(file_content[year_id[0]], /extract) year = string(fix(year_str[2]),format='(i4)') month_str = strsplit(file_content[month_id[0]], /extract) month = string(fix(month_str[2]),format='(i2.2)') day_str = strsplit(file_content[day_id[0]], /extract) day = string(fix(day_str[2]),format='(i2.2)') hour_str = strsplit(file_content[hour_id[0]], /extract) hour = string(fix(hour_str[2]),format='(i2.2)') minute_str = strsplit(file_content[minute_id[0]], /extract) minute = string(fix(minute_str[2]),format='(i2.2)') second_str = strsplit(file_content[second_id[0]], /extract) second = string(round(float(second_str[2])),format='(i2.2)') ; --- title (cic085_2018-11-30T08:16:00_g01) title = leo+'_'+year+'-'+month+'-'+day+'T'+hour+':'+minute+':'+second+'_'+gps+'_'+agency+'_'+mission ; --- lat/lon lat_id=where(strpos(file_content, '005001') eq 0) ; lat lon_id=where(strpos(file_content, '006001') eq 0) ; lon ; lat/lon is in title, and in each bending_angle level ; here only use lat/lon in title lat_str = strsplit(file_content[lat_id[0]], /extract) latitude = float(lat_str[2]) lon_str = strsplit(file_content[lon_id[0]], /extract) longitude = float(lon_str[2]) ; --- bending_angle/impact_parameter roc_str = strsplit(file_content[ELRC_id[0]], /extract) roc = float(roc_str[2])/1000. impact_id=where(strpos(file_content, '002121 MEFR 0.0 HZ Mean frequency') eq 0, impact_num) impact_parameter=[] bending_angle=[] for impact_loop = 0, impact_num-1 do begin impact_parameter_str=strsplit(file_content[impact_id[impact_loop]+1], /extract) ;M IMPACT PARAMETER bending_angle_str=strsplit(file_content[impact_id[impact_loop]+2], /extract) ; RAD BENDING ANGLE if impact_parameter_str[2] eq 'MISSING' or bending_angle_str[2] eq 'MISSING' then continue impact_parameter = [impact_parameter, float(impact_parameter_str[2])] bending_angle = [bending_angle, float(bending_angle_str[2])] endfor if n_elements(bending_angle) eq 0 then begin error_string='No bending_angle in '+input_file return_code = 1B goto, jump_end endif ; --- refractivity/height height_id=where(strpos(file_content, 'RPSEQ003 REPLICATION') ge 0, height_num) height=[] refractivity=[] for height_loop = 0, height_num-1 do begin height_str=strsplit(file_content[height_id[height_loop]+1], /extract) ; M HEIGHT refractivity_str=strsplit(file_content[height_id[height_loop]+2], /extract) ; N-UNITS ATMOSPHERIC REFRACTIVITY height = [height, (height_str[2] eq 'MISSING')?!values.F_nan:(float(height_str[2])/1000.)] refractivity = [refractivity, (refractivity_str[2] eq 'MISSING')?!values.F_nan:float(refractivity_str[2])] endfor if height_num eq 0 then begin height = !values.F_nan refractivity = !values.F_nan endif ; --- temperature/watervapor/pressure/geopotential geopotential_id=where(strpos(file_content, 'RPSEQ004 REPLICATION') ge 0, geopotential_num) geopotential=[] pressure=[] temperature=[] watervapor=[] for geopotential_loop = 0, geopotential_num-1 do begin geopotential_str=strsplit(file_content[geopotential_id[geopotential_loop]+1], /extract) ; GPM GEOPOTENTIAL HEIGHT pressure_str=strsplit(file_content[geopotential_id[geopotential_loop]+2], /extract) ;PASCALS PRESSURE temperature_str=strsplit(file_content[geopotential_id[geopotential_loop]+3], /extract) ; KELVIN TEMPERATURE watervapor_str=strsplit(file_content[geopotential_id[geopotential_loop]+4], /extract) ; KG/KG SPECIFIC HUMIDITY if temperature_str[2] eq 'MISSING' or geopotential_str[2] eq 'MISSING' then continue ; height_geop = [height_geop, (geopotential_str[2] eq 'MISSING')?!values.F_nan:$ ; float((1587301.588d - sqrt(2.51952633d12 - double(geopotential_str[2])/1000.*9806.17d0/(3.079137d-6*(1.d0-0.00259d0*cos(2.d0*!pi/180.0*latitude)))))/1000.)] geopotential = [geopotential, (geopotential_str[2] eq 'MISSING')?!values.F_nan:$ float(geopotential_str[2])/1000.] pressure = [pressure, (pressure_str[2] eq 'MISSING')?!values.F_nan:(float(pressure_str[2])/100.)] ; pa -> hpa temperature = [temperature, (temperature_str[2] eq 'MISSING')?!values.F_nan:float(temperature_str[2])] if watervapor_str[2] eq 'MISSING' or pressure_str[2] eq 'MISSING' then begin watervapor = [watervapor, !values.F_nan] endif else begin ; SPECIFIC HUMIDITY (KG/KG) --> water vapor pressure (hPa) watervapor = [watervapor, float(double(pressure_str[2])/100.*double(watervapor_str[2])/(0.622+0.378*double(watervapor_str[2])))] endelse endfor if n_elements(temperature) eq 0 then begin geopotential = !values.F_nan pressure = !values.F_nan temperature = !values.F_nan watervapor = !values.F_nan endif height_geop = geopotential2geometric(latitude, geopotential) ; ---- vertical height_ba_structure={data: impact_parameter/1000. - roc, $ name: 'Impact Height, km'} height_structure={data: height, $ name: 'Height, km'} height_geop_structure={data: height_geop, $ name: 'Height, km'} ; ---- variable bending_angle_structure={data: bending_angle, $ name: 'Bending Angle, rad', $ vertical: 'Height_ba'} refractivity_structure={data: refractivity, $ name: 'Refractivity, N-units', $ vertical: 'Height'} pressure_structure={data: pressure, $ name: 'Pressure, hPa', $ vertical: 'Height_geop'} temperature_structure={data: temperature, $ name: 'Temperature, K', $ vertical: 'Height_geop'} water_vapor_structure={data: watervapor, $ name: 'Water vapor pressure, hPa', $ vertical: 'Height_geop'} time = julday(month, day, year, hour, minute, second) ; --- data_this_profile = {title: title, $ is_bad: is_bad, $ lat: latitude, $ lon: longitude, $ time: time, $ leo: leo, $ gnss: gps, $ irs: irs, $ ; --- Height_ba: ptr_new(height_ba_structure), $ Height: ptr_new(height_structure), $ Height_geop: ptr_new(height_geop_structure), $ ; --- Temperature: ptr_new(temperature_structure), $ WaterVapor: ptr_new(water_vapor_structure), $ Pressure: ptr_new(pressure_structure), $ Refractivity: ptr_new(refractivity_structure), $ BendingAngle: ptr_new(bending_angle_structure) $ } data_this_file =[data_this_file, data_this_profile] return_code = 2B jump_end: if return_code eq 0B then begin cmd=string('echo -e "'+error_string+'" | mail -s "Alert: something is wrong in Extract_UCARbufr.pro " xinjia.zhou@noaa.gov') spawn, cmd endif RETURN, return_code END