; ---------------------------------------------------------------------------------------------------- ; This code extracts/plots the timeseries of TEMPO true and derived tropospheric columns and ; WRF VMR and emissions of NO2 ; ;- Procedures Used - ; read_TEMPO_hourly.pro ; plot_scatter_tempo_wrf.pro ; This plots the scatter plot ; plot_trop_col_wrf.pro ; This plots the spatial distribution ; con_lc.pro ; ; In this code there are multiple layers ; (1) Species ; (2) Location [lat/lon], grids ; (3) Difference / Single time-series ; (3A) Single -> 1 control scenario ; Time series of Derived/True column/Total Emissions (NO2 only) ; (3B) Difference -> 1 control/difference scenario ; Time series of Del[=(Diff. SCN - control SCN)] columns ; ; Please note: emission timeseries have same dates as TEMPO data (12-23 hours each day) ; ---------------------------------------------------------------------------------------------------- @read_TEMPO_hourly.pro !quiet=1 ; ------------------------- Paths and relevant Info --------------------------- species = ['GLYX','H2O','HCHO','NO2','O3','SO2'] tmp_path = '/home/pfister/TEMPO/FRAPPE_TEMPO/' wrf_path = '/home/pfister/TEMPO/FRAPPE_WRFchem/TEMPO/' out_path = '/home/pfister/TEMPO/Piyush/ts_out/' wrf_stup = ['WRFV4_S05_d02_MOZAIC_final_adj_EPA2014_newWRF_2xOG/', $ 'WRFV4_S05_d02_MOZAIC_final_adj_EPA2014_newWRF_phys/', $ 'WRFV4_S05_d02_MOZAIC_final_adj_EPA2014_newWRF_RAQMS/',$ 'WRFV4_S05_d02_MOZAIC_final_NEI2011adj_newWRF/', $ 'WRFV4_S05_d02_MOZAIC_final_S05_newWRF/' ] tem_prfx = 'l2_synth/TEMPO_201408' wrf_prfx = 'wrfout_d02_2014-08-' scenario = ['adj_EPA2014_newWRF_2xOG','adj_EPA2014_newWRF_phys', $ 'adj_EPA2014_newWRF_RAQMS','NEI2011adj_newWRF','S05_newWRF'] ; ----------------------------------------------------------------------------- ;goto, skip1 ; (1) Specify species, control scenarios, difference runs,for columns; Read input from the terminal: spc = '' READ, spc, PROMPT='Enter one species of interest e.g., [GLYX, H2O, HCHO, NO2, O3, SO2]: ' print,' ' ; (2) Inputs for lat/lon/number of grid points print,' Now please provide the long/lat/number of grid points ' b_lon = '' READ, b_lon, PROMPT=' longitude ' print,' ' b_lat = '' READ, b_lat, PROMPT=' latitude ' print,' ' nb = '' READ, nb, PROMPT=' number of bounding grid points [(could be any number starting from'+$ ' 0 (means gridcell [ii,jj] value) '+$ ' or >= 1 (means a bounding box with one or more grid pointis at each side)]' print,' ' ; (3) Enter control run and specify the differences scenario (if needed) diff = '' READ, diff, PROMPT='Do you need to calculate time series differences or extract single time series; '$ +' if "diff" is chosen user need to provide control and difference scenarios '$ +' if "only" is chosen only the control scenario is needed; Choose one "diff" or "only" ' print,' ' ; (3B/A) cont_scn = '' READ, cont_scn, PROMPT='Enter control scenario simulation prefix (5 options); '$ +'Enter 1-5 for "adj_EPA2014_newWRF_2xOG", "adj_EPA2014_newWRF_phys", '$ +'"adj_EPA2014_newWRF_RAQMS", "NEI2011adj_newWRF", "S05_newWRF" ' print,' ' ; (3B) if diff eq 'diff' then begin diff_scn = '' READ, diff_scn, PROMPT='Difference Scenarios you can choose one of the scenarios'$ +' (other than control SCN) Enter 1-5 for "adj_EPA2014_newWRF_2xOG",'$ +' "adj_EPA2014_newWRF_phys", "adj_EPA2014_newWRF_RAQMS", "NEI2011adj_newWRF",'$ +' "S05_newWRF" or use "all" for all at once' print,' ' endif ; if difference time series is chosen skip1: ; Used the section below for testing -- ; spc = 'NO2' ; b_lat = '38.5' ; b_lon = '-105.5' ; nb = '0' ; diff = 'diff' ;cont_scn = '5' ;diff_scn = '1' index = where(species eq spc) n = long(cont_scn)-1 if diff eq 'diff' then m = long(diff_scn)-1 ; ====================================================================================================== ; Finding number of TEMPO files for selected control scenario and species file_ct = strcompress(tmp_path+wrf_stup[n]+tem_prfx+'*h_CO_'+species[uint(index[0])] $ +'.nc',/remove_all) files1 = file_search(file_ct, count=numfiles1) dd1 = strmid(files1, strlen(tmp_path+wrf_stup[n]+tem_prfx), 5) derv_col_ts_c = fltarr(numfiles1) true_col_ts_c = fltarr(numfiles1) emis_tot_ts_c = fltarr(numfiles1) ; Finding number of TEMPO files for selected difference scenario and species if diff eq 'diff' then begin file_df = strcompress(tmp_path+wrf_stup[m]+tem_prfx+'*h_CO_'+species[uint(index[0])] $ +'.nc',/remove_all) files2 = file_search(file_df, count=numfiles2) dd2 = strmid(files2, strlen(tmp_path+wrf_stup[m]+tem_prfx), 5) derv_col_ts_d = fltarr(numfiles2) true_col_ts_d = fltarr(numfiles2) endif b_lat1 = float(b_lat) b_lon1 = float(b_lon) ; ================================================================================================== ; Part-3B/A : Control scenario only ; -------------------------------------------------------------------------------------------------- print, 'Reading TEMPO control scenario' ;- Reading TEMPO control scenario for nf = 0, numfiles1-1 do begin ;for nf = 0, 0 do begin filename = strcompress(tmp_path+wrf_stup[n]+tem_prfx+dd1[nf]+'h_CO_'+species[uint(index[0])] $ +'.nc',/remove_all) read_TEMPO_hourly, filename, output ;print,nf,' Cont ',strmid(filename,strlen(tmp_path+wrf_stup[n]+tem_prfx),20) lat_bnd1 = output.lat_bnd lon_bnd1 = output.lon_bnd date_all1 = output.year+output.month+output.day time_all1 = output.hour+ ':'+ output.minute+ ':'+output.seconds lat_all1 = output.lat lon_all1 = output.lon trop_column_all1 = output.trop_column true_trop_column_all1 = output.true_trop_column trop_column_all1[where(trop_column_all1 lt -1.0e29)] = !Values.F_NAN true_trop_column_all1[where(true_trop_column_all1 lt -1.0e29)] = !Values.F_NAN dims = size(lat_all1) xtrack = dims[1] mstep = dims[2] ;- Reading emissions file if spc eq 'NO2' then begin ; NO2 emissions only filename_em = strcompress(wrf_path+wrf_stup[n]+wrf_prfx+dd1[nf]+':00:00',/remove_all) id = ncdf_open(filename_em,/nowrite) ncdf_varget, id, 'E_NO2', emi_no2 ncdf_varget, id, 'XLONG', wrf_lon ncdf_varget, id, 'XLAT', wrf_lat ncdf_close, id wlon = reform(wrf_lon[*,0]) nwlon = n_elements(wlon) wlat = reform(wrf_lat[0,*]) nwlat = n_elements(wlat) emi_tot = fltarr(nwlon, nwlat) ;print,nf,' Emis ',strmid(filename_em,strlen(wrf_path+wrf_stup[n]+wrf_prfx),20) ;- getting total emissions for ii = 0, nwlon-1 do begin ; lon for jj = 0, nwlat-1 do begin ; lat emi_tot[ii,jj] = total(reform(emi_no2[ii,jj,*]),/NAN) endfor ; lat endfor ; lon endif ; NO2 emissions only ; ----------------------------------------------------------------------------- ; Column/Emissions value at [ii,jj] nb = long(nb) if nb eq 0 then begin ; nb (nb=0) ; Column for ii = 0, xtrack-1 do begin ; lon for jj = 0, mstep-1 do begin ; lat mn_lat = min(reform(lat_bnd1[*,ii,jj]) ,/NAN) mx_lat = max(reform(lat_bnd1[*,ii,jj]) ,/NAN) mn_lon = min(reform(lon_bnd1[*,ii,jj]) ,/NAN) mx_lon = max(reform(lon_bnd1[*,ii,jj]) ,/NAN) ind = where((b_lon1 lt mx_lon) and (b_lon1 ge mn_lon) and $ (b_lat1 lt mx_lat) and (b_lat1 ge mn_lat), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',mn_lat,' <',b_lat,' <',mx_lat,' ',mn_lon,' <',b_lon,' <',mx_lon lon_index = ii lat_index = jj endif endfor ; lat endfor ; lon derv_col_ts_c[nf] = trop_column_all1[lon_index, lat_index] true_col_ts_c[nf] = true_trop_column_all1[lon_index, lat_index] ; Emissions (NO2 only) if spc eq 'NO2' then begin ; NO2 only for ii = 0, nwlon-2 do begin ; lon for jj = 0, nwlat-2 do begin ; lat ind = where((b_lon1 lt wlon[ii+1]) and (b_lon1 ge wlon[ii]) and $ (b_lat1 lt wlat[jj+1]) and (b_lat1 ge wlat[jj]), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',wlat[jj],' <',b_lat1,' <',wlat[jj+1],' ',wlon[ii],' <',b_lon1,' <',wlon[jj+1] wlon_index = ii wlat_index = jj endif endfor ; lat endfor ; lon emis_tot_ts_c[nf] = emi_tot[wlon_index, wlat_index] endif ; NO2 only endif ; nb (nb=0) ; ----------------------------------------------------------------------------- ; Column/Emissions averaged over a Box (if chosen) Avg.[ii-nb:ii+nb ; jj-nb:jj+nb] nb = long(nb) if nb gt 0 then begin ; nb(nb>0) ; Column for ii = 0, xtrack-1 do begin ; lon for jj = 0, mstep-1 do begin ; lat mn_lat = min(reform(lat_bnd1[*,ii,jj]) ,/NAN) mx_lat = max(reform(lat_bnd1[*,ii,jj]) ,/NAN) mn_lon = min(reform(lon_bnd1[*,ii,jj]) ,/NAN) mx_lon = max(reform(lon_bnd1[*,ii,jj]) ,/NAN) ind = where((b_lon1 lt mx_lon) and (b_lon1 ge mn_lon) and $ (b_lat1 lt mx_lat) and (b_lat1 ge mn_lat), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',mn_lat,' <',b_lat,' <',mx_lat,' ',mn_lon,' <',b_lon,' <',mx_lon lon_index = ii lat_index = jj endif endfor ; lat endfor ; lon derv_col_ts_c[nf] = mean(reform(trop_column_all1[lon_index-nb:lon_index+nb, lat_index-nb:lat_index+nb]),/NAN) true_col_ts_c[nf] = mean(reform(true_trop_column_all1[lon_index-nb:lon_index+nb, lat_index-nb:lat_index+nb]),/NAN) ; Emissions if spc eq 'NO2' then begin ; NO2 only for ii = 0, nwlon-2 do begin ; lon for jj = 0, nwlat-2 do begin ; lat ind = where((b_lon1 lt wlon[ii+1]) and (b_lon1 ge wlon[ii]) and $ (b_lat1 lt wlat[jj+1]) and (b_lat1 ge wlat[jj]), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',wlat[jj],' <',b_lat1,' <',wlat[jj+1],' ',wlon[ii],' <',b_lon1,' <',wlon[jj+1] wlon_index = ii wlat_index = jj endif endfor ; lat endfor ; lon emis_tot_ts_c[nf] = mean(reform(emi_tot[wlon_index-nb:wlon_index+nb, wlat_index-nb:wlat_index+nb]), /NAN) endif ; NO2 only endif ; nb (nb>0) endfor ; numfiles1 ; ================================================================================================== ; Part-3B : If difference is chosen, this section of the code also works ; -------------------------------------------------------------------------------------------------- print, 'Reading TEMPO Difference scenario' if diff eq 'diff' then begin ; Reading TEMPO difference scenario for nf = 0, numfiles2-1 do begin ; numfiles2 filename = strcompress(tmp_path+wrf_stup[m]+tem_prfx+dd1[nf]+'h_CO_'+species[uint(index[0])] $ +'.nc',/remove_all) read_TEMPO_hourly, filename, output ;print,nf,' Diff ',strmid(filename,strlen(tmp_path+wrf_stup[m]+tem_prfx),20) lat_bnd1 = output.lat_bnd lon_bnd1 = output.lon_bnd date_all1 = output.year+output.month+output.day time_all1 = output.hour+ ':'+ output.minute+ ':'+output.seconds lat_all1 = output.lat lon_all1 = output.lon trop_column_all1 = output.trop_column true_trop_column_all1 = output.true_trop_column trop_column_all1[where(trop_column_all1 lt -1.0e29)] = !Values.F_NAN true_trop_column_all1[where(true_trop_column_all1 lt -1.0e29)] = !Values.F_NAN dims = size(lat_all1) xtrack = dims[1] mstep = dims[2] ; ----------------------------------------------------------------------------- ; Column value at [ii,jj] nb = long(nb) if nb eq 0 then begin ; nb (nb=0) ; Column for ii = 0, xtrack-1 do begin ; lon for jj = 0, mstep-1 do begin ; lat mn_lat = min(reform(lat_bnd1[*,ii,jj]) ,/NAN) mx_lat = max(reform(lat_bnd1[*,ii,jj]) ,/NAN) mn_lon = min(reform(lon_bnd1[*,ii,jj]) ,/NAN) mx_lon = max(reform(lon_bnd1[*,ii,jj]) ,/NAN) ind = where((b_lon1 lt mx_lon) and (b_lon1 ge mn_lon) and $ (b_lat1 lt mx_lat) and (b_lat1 ge mn_lat), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',mn_lat,' <',b_lat,' <',mx_lat,' ',mn_lon,' <',b_lon,' <',mx_lon lon_index = ii lat_index = jj endif endfor ; lat endfor ; lon derv_col_ts_d[nf] = trop_column_all1[lon_index, lat_index] true_col_ts_d[nf] = true_trop_column_all1[lon_index, lat_index] endif ; nb (nb = 0) ; ----------------------------------------------------------------------------- ; Columns averaged over a Box (if chosen) Avg.[ii-nb:ii+nb ; jj-nb:jj+nb] nb = long(nb) if nb gt 0 then begin ; nb (nb>0) ; Column for ii = 0, xtrack-1 do begin ; lon for jj = 0, mstep-1 do begin ; lat mn_lat = min(reform(lat_bnd1[*,ii,jj]) ,/NAN) mx_lat = max(reform(lat_bnd1[*,ii,jj]) ,/NAN) mn_lon = min(reform(lon_bnd1[*,ii,jj]) ,/NAN) mx_lon = max(reform(lon_bnd1[*,ii,jj]) ,/NAN) ind = where((b_lon1 lt mx_lon) and (b_lon1 ge mn_lon) and $ (b_lat1 lt mx_lat) and (b_lat1 ge mn_lat), nele) if ind[0] ne -1 then begin ;print, ii,' ',jj,' YES ',mn_lat,' <',b_lat,' <',mx_lat,' ',mn_lon,' <',b_lon,' <',mx_lon lon_index = ii lat_index = jj endif endfor ; lat endfor ; lon derv_col_ts_d[nf] = mean(reform(trop_column_all1[lon_index-nb:lon_index+nb, lat_index-nb:lat_index+nb]),/NAN) true_col_ts_d[nf] = mean(reform(true_trop_column_all1[lon_index-nb:lon_index+nb, lat_index-nb:lat_index+nb]),/NAN) endif ; nb (nb>0) endfor ; numfiles2 ;------ Getting the difference time series ------- derv_col_ts_df = derv_col_ts_d - derv_col_ts_c true_col_ts_df = true_col_ts_d - true_col_ts_c ; ------------------------------------------------ endif ; diff scenario ; =================================================================================================== ; writing output files (time series) 3 for single scenario and 1 for diffrence scenario ; --------------------------------------------------------------------------------------------------- print, 'Writting oputput time series files ' ; --------------------------------------------------------- ; Column file ; --------------------------------------------------------- if diff eq 'only' then begin outfile = strcompress(out_path+'ts_'+species[uint(index[0])]+'_Colm_'+scenario[n]+'_'+$ b_lon+'_'+b_lat+'_'+string(nb)+'.nc' ,/remove_all) ;spawn, 'rm -rf ' + outfile id = ncdf_create(outfile ,/clobber) xdim = ncdf_dimdef(id, 'n_hours', numfiles1) tru_id = ncdf_vardef(id, 'true_column', [xdim], /float) der_id = ncdf_vardef(id, 'derv_column', [xdim], /float) ncdf_control, id, /endef ncdf_varput, id, der_id, derv_col_ts_c ncdf_varput, id, tru_id, true_col_ts_c ncdf_control, id, /redef ncdf_attput, id, /global, "Title", "TEMPO_COLUMNs_Derived&True" ncdf_attput, id, /global, "Missing_Values", "NaNf" ncdf_attput, id, /global, "Column_units", "mol/cm2" ncdf_close, id endif ; --------------------------------------------------------- ; Emisisons file (NO2 only) ; --------------------------------------------------------- if spc eq 'NO2' then begin ; NO2 emissions only outfile = strcompress(out_path+'ts_'+species[uint(index[0])]+'_Emis_'+scenario[n]+'_'+$ b_lon+'_'+b_lat+'_'+string(nb)+'.nc' ,/remove_all) ;spawn, 'rm -rf ' + outfile id = ncdf_create(outfile ,/clobber) xdim = ncdf_dimdef(id, 'n_hours', numfiles1) emi_id = ncdf_vardef(id, 'emis_total', [xdim], /float) ncdf_control, id, /endef ncdf_varput, id, emi_id, emis_tot_ts_c ncdf_control, id, /redef ncdf_attput, id, /global, "Title", "NO2_Emissions_total" ncdf_attput, id, /global, "Missing_Values", "NaNf" ncdf_attput, id, /global, "NO2_Emissions_units", "mol/km2/hr" ncdf_close, id endif ; --------------------------------------------------------- ; Difference file ; --------------------------------------------------------- if diff eq 'diff' then begin ;spawn, 'rm -rf ' + outfile outfile = strcompress(out_path+'ts_'+species[uint(index[0])]+'_Colm_'+scenario[n]+'-'$ +scenario[m]+'_'+b_lon+'_'+b_lat+'_'+string(nb)+'.nc' ,/remove_all) id = ncdf_create(outfile ,/clobber) xdim = ncdf_dimdef(id, 'n_hours', numfiles1) tru_id = ncdf_vardef(id, 'true_column', [xdim], /float) der_id = ncdf_vardef(id, 'derv_column', [xdim], /float) ncdf_control, id, /endef ncdf_varput, id, der_id, derv_col_ts_df ncdf_varput, id, tru_id, true_col_ts_df ncdf_control, id, /redef ncdf_attput, id, /global, "Title", "TEMPO_COLUMNs_Derived&True" ncdf_attput, id, /global, "Missing_Values", "NaNf" ncdf_attput, id, /global, "Column_units", "mol/cm2" ncdf_close, id endif ; --------------------------------------------------------------------------------------------------- END