FUNCTION SUB_gpsro_qsum_ba_daily, thisDATE, mission_this, model_list ; call by gpsro_qsum_ba.pro print,thisDATE, ' SUB_gpsro_qsum_ba_daily in '+thisDATE+' begin at ',systime() ; --- folder-name/file-name use publisher (GeoOptics/S4bufr) file_pattern = 'input/'+mission_this+'/*_noref/*_'+thisDATE+'.sav' ;file_pattern = 'input/'+mission_this+'/UCARwet_noref/*_'+thisDATE+'.sav' ;file_pattern = 'input/'+mission_this+'/IGOR_noref/*_'+thisDATE+'.sav' file_list=file_search(file_pattern, count=cnt) if cnt eq 0 then begin print,'there is no files at: '+file_pattern return, 0 endif ; ------- ; initial ; is_new_want = 0 ; manual, for using the existing .sav file is_new_want = 1 ; default, generate new .sav file filename='./input_STAT/'+mission_this+'/qsum_'+thisDATE+'_ba.sav' if file_test(filename) and is_new_want eq 0 then begin restore, filename ; --- read data_prof --- variable_list = ['BA','Refractivity'] for variable_loop = 0, n_elements(variable_list)-1 do begin variable = variable_list[variable_loop] for model_loop = 0, n_elements(model_list)-1 do begin model = model_list[model_loop] ; --- in height status = execute(variable+'_total_h_'+model+' = *(data_prof_h_'+model+'.'+variable+'_total)') status = execute(variable+'_square_h_'+model+' = *(data_prof_h_'+model+'.'+variable+'_square)') status = execute(variable+'_num_h_'+model+' = *(data_prof_h_'+model+'.'+variable+'_num)') status = execute(variable+'_combine_h_'+model+' = *(data_prof_h_'+model+'.'+variable+'_combine)') endfor endfor ; variable_loop is_new_result=0 endif else begin ; is_new_want = 1 or filename non-exist is_new_result=1 data_prof={ $ ; --------------------------------------------------- ba_combine: ptr_new(/ALLOCATE_HEAP),$ ; mission(publisher) ba_num: ptr_new(/ALLOCATE_HEAP),$ ba_total: ptr_new(/ALLOCATE_HEAP),$ ba_square: ptr_new(/ALLOCATE_HEAP),$ ; --------------------------------------------------- Refractivity_combine: ptr_new(/ALLOCATE_HEAP),$ Refractivity_num: ptr_new(/ALLOCATE_HEAP),$ Refractivity_total: ptr_new(/ALLOCATE_HEAP),$ Refractivity_square: ptr_new(/ALLOCATE_HEAP) } data_prof_h_gfs = data_prof data_prof_h_gdas = data_prof data_prof_h_era = data_prof data_prof_h_era5f = data_prof endelse for model_loop = 0, n_elements(model_list)-1 do begin if model_list[model_loop] eq 'gfs' then begin status = SUB_INPUT_GFS( thisDATE, $ ; ------- lon=lon_gfs, $ ; [ 0 ~ 360] lat=lat_gfs, $ ; [90 ~ -90] ; ------- gfs_geop_vector, $ gfs_t_vector, $ ; 31:5:360:181, height:time:lon:lat gfs_wv_vector) endif else if model_list[model_loop] eq 'gdas' then begin status = SUB_INPUT_GFS( thisDATE, $ ; ------- lon=lon_gdas, $ ; [ 0 ~ 360] lat=lat_gdas, $ ; [90 ~ -90] ; ------- gdas_geop_vector, $ gdas_t_vector, $ ; 31:5:360:181, height:time:lon:lat gdas_wv_vector, $ type = 'GDAS') endif else if model_list[model_loop] eq 'era' then begin status = SUB_INPUT_ERA5_cds( thisDATE, $ ; ------- lon=lon_era, $ ; [ 0 ~ 360] lat=lat_era, $ ; [90 ~ -90] ; ------- era_geop_vector, $ era_t_vector, $ ; era5: 37:5:1441:721, height:time:lon:lat era_wv_vector) endif else if model_list[model_loop] eq 'era5f' then begin status = SUB_INPUT_ERA5_forecast( thisDATE, $ ; ------- lon=lon_era5f, $ ; [ 0 ~ 360] lat=lat_era5f, $ ; [90 ~ -90] ; ------- era5f_geop_vector, $ era5f_t_vector, $ ; era5: 37:5:1441:721, height:time:lon:lat era5f_wv_vector) endif endfor ; model_list ;for file_loop = 0, cnt-1 do begin for file_loop = cnt-1, 0, -1 do begin ; make sure UCAR is the last one ;for file_loop = 0, 0 do begin ; --- there is no fm file for S4bufr --- if strpos(file_list[file_loop], 'S4bufr') ge 0 then continue ; --- in bufr file, BA and HEIGHT are not match --- ;if strpos(file_list[file_loop], 'bufr') ge 0 then begin ; ;print, file_list[file_loop] ; is_bufr = 1 ; ;continue ;endif else is_bufr = 0 ; publisher ;file_list=file_search('input/'+mission_this+'/*_noref/*_'+thisDATE+'.sav', count=cnt) dirname = strsplit(file_dirname(file_list[file_loop]), '/', /extract) publisher = strmid(dirname[-1], 0, strlen(dirname[-1])-6) ; UCARdry/ROPPdry etc for individual RO basename = strsplit(file_basename(file_list[file_loop]), '_', /extract) satellite = basename[0] satellite2 = satellite if strpos(mission_this, 'STAR_') eq 0 or strpos(mission_this, 'UCAR_atmPrf') eq 0 then begin mission_fm = publisher endif else begin mission_array = strsplit(mission_this, '_', /extract) mission_fm = mission_array[0] if strpos(mission_this, '_nrt') gt 0 then mission_fm=mission_fm+'_nrt' $ else if strpos(mission_this, '_v2021') gt 0 then mission_fm=mission_fm+'_v2021' endelse mission_fm2 = mission_fm if strpos(mission_this, 'ROMEX_STAR_ROPP') eq 0 or strpos(mission_this, 'ROMEX_STAR_RFSI') eq 0 then begin if publisher eq 'GEOOPTICS' then begin mission_fm = 'GeoOptics_nrt' satellite2 = 'geooptics' endif else if publisher eq 'IGOR' and satellite2 eq 'tsx' then begin satellite = 'terrasarx' endif else if publisher eq 'IGOR' and satellite2 eq 'tdx' then begin satellite = 'tandemx' endif else if publisher eq 'METOP' then begin mission_fm = 'GRAS' satellite2 = 'metop' endif else begin mission_fm = publisher endelse endif combine_flag = 0 if mission_this eq 'COSMIC' or strpos(mission_this, 'COSMIC_') eq 0 or $ mission_this eq 'COSMIC1' or strpos(mission_this, 'COSMIC1_') eq 0 then begin combine_this='COSMIC1 ('+publisher+')' endif else if mission_this eq 'COSMIC2' or strpos(mission_this, 'COSMIC2_') eq 0 then begin combine_this='COSMIC2 ('+publisher+')' endif else if strpos(mission_this, 'SPIRE') eq 0 then begin combine_this='Spire ('+publisher+')' endif else if strpos(mission_this, 'PlanetiQ') eq 0 then begin combine_this='PlanetiQ ('+publisher+')' endif else if strpos(mission_this, 'GeoOptics') eq 0 then begin combine_this='GeoOptics ('+publisher+')' endif else if strpos(mission_this, 'UCAR_') eq 0 or $ strpos(mission_this, 'STAR_') eq 0 or $ strpos(mission_this, 'ROMEX_') eq 0 then begin if strmid(publisher, 3, /REVERSE_OFFSET) eq '_nrt' then begin combine_this=satellite+'rt' endif else begin if satellite eq 'tsx' then combine_this = 'terrasarx' $ else if satellite eq 'tdx' then combine_this = 'tandemx' $ else if satellite eq 'cicero' then combine_this = 'geooptics' $ else if strpos(mission_this, 'ROMEX_') eq 0 and satellite eq 'metopa' then combine_this = 'metop' $ else if strpos(mission_this, 'ROMEX_') eq 0 and satellite eq 'metopb' then combine_this = 'metop' $ else if strpos(mission_this, 'ROMEX_') eq 0 and satellite eq 'metopc' then combine_this = 'metop' $ else combine_this = satellite endelse endif else begin ;combine_this=leo+' ('+publisher+')' ; separate each leo satellite, like MetopA/MetopB, or TerrasaX/TendenX combine_flag = 1 endelse ; --------------------- file_gfs ='input_fm/'+mission_fm+'/UCARdry_STAR_gfs06/'+satellite+'_'+thisDATE+'.sav' file2_gfs ='input_fm/'+mission_fm2+'/UCARdry_STAR_gfs06/'+satellite2+'_'+thisDATE+'.sav' file_gdas ='input_fm/'+mission_fm+'/UCARdry_STAR_gfs00/'+satellite+'_'+thisDATE+'.sav' file2_gdas ='input_fm/'+mission_fm2+'/UCARdry_STAR_gfs00/'+satellite2+'_'+thisDATE+'.sav' file_era ='input_fm/'+mission_fm+'/UCARdry_STAR_era5/'+satellite+'_'+thisDATE+'.sav' file2_era ='input_fm/'+mission_fm2+'/UCARdry_STAR_era5/'+satellite2+'_'+thisDATE+'.sav' if strpos(publisher, '_UMD_') gt 0 then begin file_era5f='input_fm/'+mission_fm+'/ROPPdry_UMD_STAR_era5f/'+satellite+'_'+thisDATE+'.sav' file2_era5f='input_fm/'+mission_fm2+'/ROPPdry_UMD_STAR_era5f/'+satellite2+'_'+thisDATE+'.sav' endif else begin file_era5f='input_fm/'+mission_fm+'/ROPPdry_STAR_era5f/'+satellite+'_'+thisDATE+'.sav' file2_era5f='input_fm/'+mission_fm2+'/ROPPdry_STAR_era5f/'+satellite2+'_'+thisDATE+'.sav' endelse print, string(10B)+'----' print, 'publisher:', publisher print,'input_file: '+file_list[file_loop]+string(10B) if file_test(file_gfs) then begin print,'file_gfs: '+file_gfs restore,file_gfs data_gfs = data endif else if file_test(file2_gfs) then begin print,'file2_gfs: '+file2_gfs restore,file2_gfs data_gfs = data endif else begin print, 'file_gfs non-exist:'+file_gfs print, 'file2_gfs non-exist:'+file2_gfs data_gfs=[] endelse if file_test(file_gdas) then begin print,'file_gdas: '+file_gdas restore,file_gdas data_gdas = data endif else if file_test(file2_gdas) then begin print,'file2_gdas: '+file2_gdas restore,file2_gdas data_gdas = data endif else begin print, 'file_gdas non-exist:'+file_gdas print, 'file2_gdas non-exist:'+file2_gdas data_gdas=[] endelse if file_test(file_era) then begin print,'file_era: '+file_era restore,file_era data_era = data endif else if file_test(file2_era) then begin print,'file2_era: '+file2_era restore,file2_era data_era = data endif else begin print, 'file_era non-exist:'+file_era print, 'file2_era non-exist:'+file2_era data_era = [] endelse if file_test(file_era5f) then begin print,'file_era5f: '+file_era5f restore,file_era5f data_era5f = data endif else if file_test(file2_era5f) then begin print,'file2_era5f: '+file2_era5f restore,file2_era5f data_era5f = data endif else begin print, 'file_era5f non-exist:'+file_era5f print, 'file2_era5f non-exist:'+file2_era5f data_era5f = [] endelse if n_elements(data_gfs) eq 0 and $ n_elements(data_gdas) eq 0 and $ n_elements(data_era) eq 0 and $ n_elements(data_era5f) eq 0 then begin print, 'All fm files non-exist' continue endif ;print, 'size(data_gfs)',size(data_gfs) ;print, 'size(data_era)',size(data_era) ;print, 'size(data_era5f)',size(data_era5f) ; --------------------- delvar, data ; only work in main level ;undefine, data ; ;abc=temporary(data) ;print,file_list[file_loop] restore,file_list[file_loop] ;print, 'size(data)',size(data) if total(tag_names(data) eq 'BENDINGANGLE') eq 0 then begin print, "total(tag_names(data) eq 'BENDINGANGLE') eq 0" continue endif if total(tag_names(data) eq 'HEIGHT') eq 0 then begin print, "total(tag_names(data) eq 'HEIGHT') eq 0" continue ; JPLdry have no HEIGHT, but HEIGHT_BA. It's doable, not require a lot of updates endif data = data[where(data.is_bad eq 0, /null)] ;print, 'size(data)',size(data) non_match = lonarr(n_elements(model_list)) for profile_loop = 0, n_elements(data)-1 do begin ;print, 'profile_loop:',profile_loop title_full = strsplit(data[profile_loop].title,'_', /extract) hour = strmid(title_full[1],11, 2) minute= strmid(title_full[1],14, 2) time_loc_2 = (hour+minute/60.)/12. ; time location for era5_forecast, 12 hours interval, 2 reports per day (06:00/18:00) time_loc_4 = (hour+minute/60.)/6. ; time location for era5/gfs/gdas, 6 hours interval, 4 reports per day time_loc_8 = (hour+minute/60.)/3. ; time location for merra, 3 hours interval, 8 reports per day leo = data[profile_loop].leo gnss = data[profile_loop].gnss lat_2 = data[profile_loop].lat lon_2 = data[profile_loop].lon time_2 = data[profile_loop].time if combine_flag eq 1 then begin combine_this=leo+' ('+publisher+')' ; separate each leo satellite, like MetopA/MetopB, or TerrasaX/TendenX endif impact_sat = *(data[profile_loop].Height_ba) impact_sat = impact_sat.data height_sat = *(data[profile_loop].Height) height_sat = height_sat.data ba_sat = *(data[profile_loop].bendingangle) ba_sat = ba_sat.data ref_sat = *(data[profile_loop].REFRACTIVITY) ref_sat = ref_sat.data if n_elements(impact_sat) ne n_elements(ba_sat) then begin print, file_list[file_loop]+' n_elements(impact_sat) ne n_elements(ba_sat)', n_elements(impact_sat), n_elements(ba_sat) continue endif ; --- interpolation --- ; for model_loop = 0, n_elements(model_list)-1 do begin ;print, 'model_loop:', model_loop model = model_list[model_loop] cmd='is_valid = n_elements(data_'+model+')' res=execute(cmd) if is_valid eq 0 then begin ;print, 'is_valid for model eq 0' continue ; model in non-exist endif res=execute('data_model = data_'+model) if model eq 'era5f' then time_loc = time_loc_2 $ else time_loc = time_loc_4 dif_time = abs(data_model.time - time_2) cos_distance = cos(lat_2/180.*!dpi)*cos(data_model.lat/180.*!dpi)*$ cos(lon_2/180.*!dpi-data_model.lon/180.*!dpi)+$ sin(lat_2/180.*!dpi)*sin(data_model.lat/180.*!dpi) cos_distance = (cos_distance <=1) distance = 6371.004*acos(cos_distance) if mission_this eq 'PlanetiQ_nrt_STAR_v2.37' then begin ; by Xinjia, the UMD leo name is 'P004'/'P005'. But in UCAR files, leo name is 'GN04'/'GN05' ro_leo = strmid(leo, 2,2) model_leo = strmid(data_model.leo, 2, 2) endif else begin ro_leo = leo model_leo = data_model.leo endelse match_id = where( dif_time lt 5./60./24. and $ ; minute -> hour -> day distance lt 200 and $ ; km strupcase(data_model.gnss) eq strupcase(gnss) and $ strupcase(model_leo) eq strupcase(ro_leo), match_num) res=execute(cmd) if match_num eq 0 then begin ;print, '(SUB_gpsro_qsum_ba_daily.pro) match_num eq 0: ',data[profile_loop].title non_match[model_loop] ++ continue ; this is possible, becase ROPP fm may stop for some reason endif ;print, data[profile_loop] ;print, data_era5f[match_id[0]] res = execute('ba_fm = *(data_'+model+'[match_id[0]].BENDINGANGLE)') ba_fm = ba_fm.data ref = execute('ref_fm = *(data_'+model+'[match_id[0]].REFRACTIVITY)') ref_fm = ref_fm.data res = execute('height_fm = *(data_'+model+'[match_id[0]].HEIGHT_LEVEL1B)') height_fm = height_fm.data ;print, 'read Height_ba from fm', model_loop res = execute('impact_fm = *(data_'+model+'[match_id[0]].HEIGHT_LEVEL1a)') impact_fm = impact_fm.data ;print, combine_this ;print, 'size(impact_sat)', size(impact_sat) ;print, 'size(ba_sat)', size(ba_sat) ;print, 'size(ref_sat)', size(ref_sat) ;print, 'size(ba_fm)', size(ba_fm) ;print, 'size(ref_fm)', size(ref_fm) ;cmd = 'print, "size(model_geop_vector)", size('+model+'_geop_vector)' ;res2=execute(cmd) ;print, 'for ba' cmd = 'SUB_interpolate_model_ba, lat_2, lon_2, lat_'+model+', lon_'+model+', time_loc, '+model+'_geop_vector, combine_this, ' +$ 'ba_sat, height_sat, impact_sat, is_ba=1, '+$ 'ba_fm, height_fm, impact_fm, '+$ ; --- 'ba_total_h_'+model+', ba_square_h_'+model+', ba_num_h_'+model+', ba_combine_h_'+model res=execute(cmd) ;print, 'for ref' cmd = 'SUB_interpolate_model_ba, lat_2, lon_2, lat_'+model+', lon_'+model+', time_loc, '+model+'_geop_vector, combine_this, ' +$ 'ref_sat, height_sat, impact_sat, is_ba=0, '+$ 'ref_fm, height_fm, impact_fm, '+$ ; --- 'Refractivity_total_h_'+model+', Refractivity_square_h_'+model+', Refractivity_num_h_'+model+', Refractivity_combine_h_'+model res=execute(cmd) ; print, cmd ;print, 'ba_combine_h_era5f:' ;print, ba_combine_h_era5f ;print, ba_total_h_era5f ;print, 'Refractivity_combine_h_era5f:' ;print, Refractivity_combine_h_era5f ;print, Refractivity_total_h_era5f endfor ; model ;print, 'after model' endfor ; profile ;print, 'after profile' ;print, 'Refractivity_combine_era:',Refractivity_combine_era print, '(SUB_gpsro_qsum_ba_daily.pro) non_match=' print, ' ',non_match endfor ; file ; ----------------- variable_list = ['BA','Refractivity'] for variable_loop = 0, n_elements(variable_list)-1 do begin variable = variable_list[variable_loop] for model_loop = 0, n_elements(model_list)-1 do begin model = model_list[model_loop] ; --- in height only --- cmd = 'data_prof_h_'+model+'.'+variable+'_total=ptr_new('+variable+'_total_h_'+model+')' status = execute(cmd) cmd = 'data_prof_h_'+model+'.'+variable+'_square=ptr_new('+variable+'_square_h_'+model+')' status = execute(cmd) cmd = 'data_prof_h_'+model+'.'+variable+'_num=ptr_new('+variable+'_num_h_'+model+')' status = execute(cmd) cmd = 'data_prof_h_'+model+'.'+variable+'_combine=ptr_new('+variable+'_combine_h_'+model+')' status = execute(cmd) endfor ; model_loop endfor ; variable_loop ;print, Temperature_combine_era_all_p ;print, Temperature_num_era_all_p file_mkdir, './input_STAT/'+mission_this save, data_prof_h_gfs, $ data_prof_h_gdas, $ data_prof_h_era, $ data_prof_h_era5f, $ filename=filename, /COMPRESS print, 'filename:',filename print,'SUB_gpsro_qsum_ba_daily.pro done for mission:'+mission_this+' in '+thisDATE+' ',systime(), string(10B) return, 2 END