pro seac4rs_transport ;dir = '/Volumes/Data_Domino/Data/' dir = '/Volumes/Data_Gator/Data/' dc8_10s = hash() dc8_was = hash() er2_10s = hash() er2_was = hash() restore,dir+'Aircraft/Missions/seac4rs/dc8/seac4rs_dc8_10s_merge.sav' nv_dc8_10s = n_elements(varnames_dc8) for i = 0, nv_dc8_10s-1 do dc8_10s[varnames_dc8[i]] = reform(seac4rs_dc8_10s.field001[i,*]) restore,dir+'Aircraft/Missions/seac4rs/dc8/seac4rs_dc8_was_merge.sav' nv_dc8_was = n_elements(varnames_dc8_was) for i = 0, nv_dc8_10s-1 do dc8_was[varnames_dc8_was[i]] = reform(seac4rs_dc8_was.field001[i,*]) restore,dir+'Aircraft/Missions/seac4rs/er2/seac4rs_er2_10s_merge.sav' nv_er2_10s = n_elements(varnames_er2) for i = 0, nv_er2_10s-1 do er2_10s[varnames_er2[i]] = reform(seac4rs_er2_10s.field001[i,*]) restore,dir+'Aircraft/Missions/seac4rs/er2/seac4rs_er2_was_merge.sav' nv_er2_was = n_elements(varnames_er2_was) for i = 0, nv_er2_was-1 do er2_was[varnames_er2_was[i]] = reform(seac4rs_er2_was.field001[i,*]) ;print,varnames_er2 er2_10s['JDAY'] += (julday(1,1,2013,0)-1.) + er2_10s['UTC']/86400. caldat,er2_10s['JDAY'],mon,day,yr,hr er2_was['JDAY'] += (julday(1,1,2013,0)-1.) + er2_10s['UTC']/86400. caldat,er2_was['JDAY'],mon,day,yr,hr ;print,mon[0],day[0],yr[0] ;print,er2_was['UTC',0:10] ; Read in MERRA2 tropopause. ncdf_get,dir+'/Reanalysis/MERRA2/tp.monmean.zm.nc',['lat','time','tpz','tpp','tpt','ctpt','ctpz'],tpzm,/quiet caldat,tpzm['time','value',404]+julday(1,1,1900,0),mon,day,yr tpzm_theta = tpzm['tpt','value'] * (1e3/tpzm['tpp','value'])^0.286 ncdf_get,dir+'/Reanalysis/MERRA2/tp-201308.nc',['lon','lat','time','tpz','tpp','tpt','ctpt','ctpz'],tpm,/quiet nym = (tpm['lat','dim_sizes'])[0] mer_tp = tpm['tpp','value'] mer_tpt = tpm['tpt','value'] mer_lat = tpm['lat','value'] lon = tpm['lon','value'] xi = where(lon ge -112 and lon le -85) tpm_theta = fltarr(nym) & tpm_theta_sd = tpm_theta & tpause = fltarr(nym) for y = 0, nym-1 do begin tpause[y] = mean(tpzm['tpp','value',y,*]) tp_th = mer_tpt[xi,y,80:-1] * (1e3/mer_tp[xi,y,80:-1])^0.286 stats = moment(tp_th,sdev=sdev) tpm_theta[y] = stats[0] tpm_theta_sd[y] = sdev endfor ncdf_get,dir+'/Reanalysis/MERRA2/tp-201309.nc',['lon','lat','time','tpz','tpp','tpt','ctpt','ctpz'],tpm,/quiet mer_tp = tpm['tpp','value'] mer_tpt = tpm['tpt','value'] tpm_theta2 = fltarr(nym) & tpm_theta_sd2 = tpm_theta & tp_m = tpm_theta for y = 0, nym-1 do begin tp_th = mer_tpt[xi,y,80:-1] * (1e3/mer_tp[xi,y,80:-1])^0.286 stats = moment(tp_th,sdev=sdev) tpm_theta2[y] = stats[0] tpm_theta_sd2[y] = sdev tp_m[y] = mean(mer_tp[xi,y,80:-1]) endfor restore,dir+'/Reanalysis/MERRA2/merra_zm.sav' nzm = n_elements(pssr) tpm_theta_avg = tpm_theta & tpm_theta_sd_avg = tpm_theta & theta_mean = fltarr(nym,nzm) for y = 0, nym-1 do begin tpm_theta_avg[y] = 0.5 * (tpm_theta[y] + tpm_theta2[y]) tpm_theta_sd_avg[y] = 0.5 * (tpm_theta_sd[y] + tpm_theta_sd2[y]) for z = 0, nzm-1 do theta_mean[y,z] = mean(theta_zm[y,z,*]) endfor ; Theta grid nzth = 21 dth = 10. theta_grid = findgen(nzth)*dth+295. nhist = 14 er2_was_theta_grid = replicate(!values.f_nan,3,nzth,nv_er2_was) & dc8_was_theta_grid = replicate(!values.f_nan,3,nzth,nv_dc8_was) er2_was_theta_grid_hist = replicate(!values.f_nan,nhist,nzth,nv_er2_was) & dc8_was_theta_grid_hist = replicate(!values.f_nan,nhist,nzth,nv_dc8_was) er2_was_theta_grid_hist_arr = er2_was_theta_grid_hist & dc8_was_theta_grid_hist_arr = dc8_was_theta_grid_hist for i = 0, nv_er2_was-1 do begin for z = 5, nzth-3 do begin zi = where(er2_was['THETA'] ge theta_grid[z]-dth/2.-2 and er2_was['THETA'] lt theta_grid[z]+dth/2.+2 and er2_was[varnames_er2_was[i]] gt 0.,nn) if nn gt 2 then begin si = sort(er2_was[varnames_er2_was[i],zi]) stats = moment(er2_was[varnames_er2_was[i],zi[si[1:-2]]],sdev=sdev) med = median(er2_was[varnames_er2_was[i],zi]) er2_was_theta_grid[*,z,i] = [med,stats[0],sdev] ; if varnames_er2_was[i] eq 'HCFC-142b_WAS' then print,z,er2_was[varnames_er2_was[i],zi[si]] endif zi = where(er2_was['THETA'] ge theta_grid[z]-dth and er2_was['THETA'] lt theta_grid[z]+dth and er2_was[varnames_er2_was[i]] gt 0.,nn) if nn gt 2 then begin si = sort(er2_was[varnames_er2_was[i],zi]) H = histogram(er2_was[varnames_er2_was[i],zi],nbins=nhist,max=er2_was[varnames_er2_was[i],zi[si[-1]]],min=er2_was[varnames_er2_was[i],zi[si[0]]],locations=locs) er2_was_theta_grid_hist[*,z,i] = H / total(H) er2_was_theta_grid_hist_arr[*,z,i] = locs endif endfor endfor for i = 0, nv_dc8_was-20 do begin for z = 0, nzth-3 do begin zi = where(dc8_was['THETA'] ge theta_grid[z]-dth/2. and dc8_was['THETA'] lt theta_grid[z]+dth/2. and dc8_was[varnames_dc8_was[i]] gt 0.,nn) if nn gt 2 then begin si = sort(dc8_was[varnames_dc8_was[i],zi]) stats = moment(dc8_was[varnames_dc8_was[i],zi[si[1:-2]]],sdev=sdev) med = median(dc8_was[varnames_dc8_was[i],zi]) dc8_was_theta_grid[*,z,i] = [med,stats[0],sdev] endif zi = where(dc8_was['THETA'] ge theta_grid[z]-dth and dc8_was['THETA'] lt theta_grid[z]+dth and dc8_was[varnames_dc8_was[i]] gt 0.,nn) if nn gt 2 then begin si = sort(dc8_was[varnames_dc8_was[i],zi]) H = histogram(dc8_was[varnames_dc8_was[i],zi],nbins=nhist,max=dc8_was[varnames_dc8_was[i],zi[si[-1]]],min=dc8_was[varnames_dc8_was[i],zi[si[0]]],locations=locs) dc8_was_theta_grid_hist[*,z,i] = H / total(H) dc8_was_theta_grid_hist_arr[*,z,i] = locs endif endfor endfor ; CI vars. ci_vars_er2 = ['O3_UAS','n-butane_WAS','i-butane_WAS','Propane_WAS','Ethyne_WAS','1-2-Dichloroethane_WAS','Ethane_WAS','Perchlorethylene_WAS','Chloroform_WAS', $ 'HFC-152a_WAS','MethylBromide_WAS','HCFC-141b_WAS','HCFC-22_WAS','HFC-134a_WAS','Halon-1211_WAS','HCFC-142b_WAS','CarbonTetrachloride_WAS','CFC-11_WAS','CFC-113_WAS','CFC-12_WAS','CFC-114_WAS'] ci_vars_dc8 = ['O3_ESRL','n-Butane_WAS','i-Butane_WAS','Propane_WAS','Ethyne_WAS','1_2-dichloroethane_WAS','Ethane_WAS','C2Cl4_WAS','CHCl3_WAS', $ 'HFC-152a_WAS','CH3Br_WAS','HCFC-141b_WAS','HCFC-22_WAS','HFC-134a_WAS','H-1211_WAS','HCFC-142b_WAS','CCl4_WAS','CFC-11_WAS','CFC-113_WAS','CFC-12_WAS','CFC-114_WAS'] ci_lifetimes = [100,8,9,14,16,50,80,100,200, $ 600,800,3e3,4e3,5e3,6e3,7e3,1e4,1.5e4,2e4,2.5e4,6e4] nciv = n_elements(ci_vars_er2) ; Atmospheric lifetimes from sources. ; n-butane ? days (OH reaction) WMO, 2018; 1.9 days tropics (Luo et al.); 6.8 days (Hodnebrog et al, 2018); 4.2 days NH subtropics (Tang et al. 2007); 2.5-4500 yrs NO3-O3 (NAS, 1991) ; i-butane 5-11 days (OH reaction) WMO, 2018; 2.1 days tropics (Luo et al.); 4.7 days NH subtropics (Tang et al. 2007) ; Propane 10-27 days (OH reaction) WMO, 2018; 4.2 days tropics (Luo et al.); 13 days (Hodnebrog et al, 2018); 9.2 days NH subtropics (Tang et al. 2007); 2.5-4500 yrs NO3-O3 (NAS, 1991) ; Ethyne 14 days (Xiao et al., 2007) ; 6 days tropics (Luo et al.); 11 days NH subtropics (Tang et al. 2007); 2-6 years NO3-O3 (NAS, 1991) ; Ethane 58 days (Hodnebrog et al, 2018); 18 days tropics (Luo et al.); 40 days NH subtropics (Tang et al. 2007); 12-4500 yrs NO#-O3 (NAS, 1991) ; 1-2-Dichloroethane (C2H4Cl2) 82 days (41-555 days) trop (WMO, 2018); 16 days tropics (Luo et al.) ; Perchloroethylene (C2Cl4) 111 days (Wuebbles et al., 2011); 105 days (Olaguer, 2002); 27 days tropics (Luo et al.); Cl reaction in strat appears to shorten lifetime (Singh et al., 1996) ; Chloroform (CHCl3) 6 months 2-3 months trop 1-2 yrs strat; 6 months (97-1145 days) trop (WMO, 2018); 48 days tropics (Luo et al.) ; HFC-152a 1.5 yrs ; 1.55 yrs OH trop, 39 yrs strat (WMO, 2018) ; Methyl Bromide (CH3Br) 2 yrs (trop OH reaction) 33 yrs strat (Saltzman et al., 2004); 1.8 yrs trop, 26 yrs strat (WMO, 2018) ; HCFC-141b 9.4 yrs 74 yrs strat; 10.7 yrs OH trop, 72.3 strat (WMO, 2018) ; HCFC-22 11.9 yrs 161 yrs strat; 13 yrs OH trop, 161 strat (WMO, 2018) ; HFC-134a 13.4 yrs 5 yrs tropics (Chelpon et al.), 14.1 yrs OH trop, 267 yrs strat (WMO, 2018) ; HCFC-142b 18.0 yrs 212 yrs strat; 19.3 yrs OH trop, 212 yrs strat (WMO, 2018) ; H-1211 12.0 yrs 17 yrs trop 33 yrs strat (Chipperfield et al., 2013); 16 yrs total, 41 yrs strat (WMO, 2018) ; CCl4 48 yrs 1230 yrs trop 50 yrs strat; 32 yrs total, 44 yrs strat (WMO, 2018) ; CFC-11 52 yrs (WMO, 2018) ; CFC-113 93 yrs (WMO, 2018) ; CFC-12 102 yrs (WMO, 2018) ; CFC-114 189 yrs (WMO, 2018) ; Read in NOAA surface data for boundary layer estimates. ---------------------------------------------------- restore,dir+'Surface_Trace_Gas/co2_surface_flask.sav' flask_data = replicate(!values.f_nan,nco2_mlo,6,nciv) & flask_data_ext = replicate(!values.f_nan,nco2_ext,6,nciv) for i = 20, 1, -1 do begin if i le 2 then begin stations = ['smo','mlo','key','lef'] starti = [70,70,69,69] for s = 0, 3 do begin data = read_ascii(dir+'Surface_Trace_Gas/'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'_'+stations[s]+'.txt',data_start=starti[s]) yr = reform(data.field01[1,*]) mon = reform(data.field01[2,*]) raw = reform(data.field01[11,*]) time_raw = yr + mon/12. + 1./24. time = time_raw[uniq(time_raw)] nn = n_elements(time) flask_mon = replicate(!values.f_nan,nn) for t = 0, nn-1 do begin ti = where(time_raw eq time[t] and raw gt 0,nti) if nti gt 0 then begin flask_mon[t] = mean(raw[ti]) ti = where(co2_time_ext ge time[t]-0.01 and co2_time_ext le time[t]+0.01) flask_data_ext[ti,s,i] = flask_mon[t] endif endfor month = fix(12 * (co2_time_ext - fix(co2_time_ext))) seas_cyc = fltarr(12) for t = 0, 11 do seas_cyc[t] = mean(flask_data_ext[where(month eq t and finite(flask_data_ext[*,s,i])),s,i]) chk = where(~finite(flask_data_ext[*,s,i]),nchk) for t = 0, nchk-1 do begin month = fix(12 * (co2_time_ext[chk[t]] - fix(co2_time_ext[chk[t]]))) flask_data_ext[chk[t],s,i] = seas_cyc[month] endfor flask_data_ext[*,4,i] = flask_data_ext[*,1,i] ; Tropics flask_data_ext[*,5,i] = flask_data_ext[*,3,i] ; Tropics endfor endif if ci_vars_dc8[i] eq 'Ethyne_WAS' then begin flask_data_ext[*,*,i] = flask_data_ext[*,*,i+2] / 6. t1 = where(co2_time_ext ge 2013) for y = 0, 5 do begin for t = t1[0], 0, -1 do flask_data_ext[t,y,i] *= (1. + (2013.-co2_time_ext[t])/30.) endfor flask_data_ext[*,*,i] *= 0.7 endif if i eq 3 or i eq 6 then begin data = read_ascii(dir+'Surface_Trace_Gas/HATS Zonal Means/HATS_Zonal_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.txt',data_start=11,delimiter=',') nn = n_elements(data.field01[1,*]) year = reform(data.field01[0,*]) mon = reform(data.field01[1,*]) nh_90 = reform(data.field01[2,*]) nh_60 = reform(data.field01[4,*]) nh_45 = reform(data.field01[6,*]) nh_20 = reform(data.field01[8,*]) sh_0 = reform(data.field01[10,*]) dates = year + mon/12. + 1./24. for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then begin flask_data_ext[t,5,i] = nh_45[chk] flask_data_ext[t,4,i] = nh_20[chk] flask_data_ext[t,3,i] = nh_60[chk] flask_data_ext[t,2,i] = nh_45[chk] flask_data_ext[t,1,i] = nh_20[chk] flask_data_ext[t,0,i] = sh_0[chk] endif endfor seas_cyc = fltarr(12,6) for t = 0, 11 do begin seas_cyc[t,0] = mean(sh_0[where(mon eq t+1 and finite(sh_0))]) seas_cyc[t,1] = mean(nh_20[where(mon eq t+1 and finite(nh_20))]) seas_cyc[t,2] = mean(nh_45[where(mon eq t+1 and finite(nh_45))]) seas_cyc[t,3] = mean(nh_60[where(mon eq t+1 and finite(nh_60))]) seas_cyc[t,4] = mean(nh_20[where(mon eq t+1 and finite(nh_20))]) seas_cyc[t,5] = mean(nh_45[where(mon eq t+1 and finite(nh_45))]) endfor for y = 0, 5 do begin gdy = where(finite(flask_data_ext[*,y,i])) for t = gdy[0]-1, 0, -1 do begin month = fix(12 * (co2_time_ext[t] - fix(co2_time_ext[t]))) flask_data_ext[t,y,i] = seas_cyc[month,y] if ci_vars_dc8[i] eq 'Ethane_WAS' then flask_data_ext[t,y,i] *= (1. + (2013.-co2_time_ext[t])/100.) ; Ethane growth rate from Helmig et al., 2016 endfor endfor endif if ci_vars_dc8[i] eq '1_2-dichloroethane_WAS' then begin ; flask_data_ext[*,*,i] = 13.0 flask_data_ext[*,3,i] = replicate(13.,nco2_ext) flask_data_ext[*,2,i] = replicate(13.,nco2_ext) flask_data_ext[*,1,i] = replicate(11.,nco2_ext) flask_data_ext[*,0,i] = replicate(8.,nco2_ext) endif if i ge 9 and i le 19 then begin data = read_ascii(dir+'Surface_Trace_Gas/HATS Zonal Means/HATS_Zonal_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.txt',data_start=11,delimiter=',') nn = n_elements(data.field01[1,*]) year = reform(data.field01[0,*]) mon = reform(data.field01[1,*]) nh_90 = reform(data.field01[2,*]) nh_60 = reform(data.field01[4,*]) nh_45 = reform(data.field01[6,*]) nh_20 = reform(data.field01[8,*]) sh_0 = reform(data.field01[10,*]) dates = year + mon/12. + 1./24. miss = where(~finite(nh_60),nmiss) if nmiss gt 0 then begin nh_60[miss] = replicate(mean(nh_60,/nan),nmiss) endif for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then begin flask_data_ext[t,5,i] = nh_45[chk] flask_data_ext[t,4,i] = nh_20[chk] flask_data_ext[t,3,i] = nh_60[chk] flask_data_ext[t,2,i] = nh_45[chk] flask_data_ext[t,1,i] = nh_20[chk] flask_data_ext[t,0,i] = sh_0[chk] endif endfor gd = where(finite(flask_data_ext[*,1,i])) mlo_yrsm = lowpass_cfc(flask_data_ext[gd,1,i], BOX=12, EDGE_PFCAST=1) deriv_yr = 12*sg_smooth(mlo_yrsm, NLEFT=12/2, NRIGHT=12/2-1, DERIV=1, DELTA=1.0, EDGE_PFCAST=1) mlo_gr = replicate(!values.f_nan,nco2_ext) mlo_gr[gd] = deriv_yr / mlo_yrsm for y = 0, 5 do begin gdy = where(finite(flask_data_ext[*,y,i])) for t = gdy[0]-1, 0, -1 do begin if finite(mlo_gr[t+1]) then mgr = mlo_gr[t+1] else mgr = mlo_gr[gd[0]] flask_data_ext[t,y,i] = (1.-mgr/12.) * flask_data_ext[t+1,y,i] endfor endfor ; if ci_vars_dc8[i] eq 'HFC-152a_WAS' then flask_data_ext[*,4:5,i] *= 1.05 ; if ci_vars_dc8[i] eq 'CH3Br_WAS' then flask_data_ext[*,*,i] *= 1.02 if ci_vars_dc8[i] eq 'HCFC-141b_WAS' then flask_data_ext[*,*,i] *= 1.03 if ci_vars_dc8[i] eq 'HFC-134a_WAS' then flask_data_ext[*,*,i] *= 1.05 if ci_vars_dc8[i] eq 'HCFC-142b_WAS' then flask_data_ext[*,*,i] *= 1.009 ;1.0075 if ci_vars_dc8[i] eq 'H-1211_WAS' then flask_data_ext[*,*,i] *= 1.015 ;1.005 if ci_vars_dc8[i] eq 'CCl4_WAS' then flask_data_ext[*,*,i] *= 0.98 if ci_vars_dc8[i] eq 'CFC-11_WAS' then flask_data_ext[*,*,i] *= 0.98 if ci_vars_dc8[i] eq 'CFC-113_WAS' then flask_data_ext[*,*,i] *= 1.015 endif if i eq 7 or i eq 8 or i eq 20 then begin if i eq 7 then ii = 37 if i eq 8 then ii = 46 if i eq 20 then ii = 85 data = read_ascii(dir+'Surface_Trace_Gas/AGAGE/SMO-medusa.mon.txt',data_start=16) nn = n_elements(data.field001[0,*]) dates = reform(data.field001[0,*]) flask = reform(data.field001[ii,*]) for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then flask_data_ext[t,0,i] = mean(flask[chk]) endfor data = read_ascii(dir+'Surface_Trace_Gas/AGAGE/RPB-medusa.mon.txt',data_start=16) nn = n_elements(data.field001[0,*]) dates = reform(data.field001[0,*]) flask = reform(data.field001[ii,*]) for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then flask_data_ext[t,1,i] = mean(flask[chk]) endfor data = read_ascii(dir+'Surface_Trace_Gas/AGAGE/THD-medusa.mon.txt',data_start=16) nn = n_elements(data.field001[0,*]) dates = reform(data.field001[0,*]) flask = reform(data.field001[ii,*]) for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then flask_data_ext[t,2,i] = mean(flask[chk]) endfor data = read_ascii(dir+'Surface_Trace_Gas/AGAGE/MHD-medusa.mon.txt',data_start=16) nn = n_elements(data.field001[0,*]) dates = reform(data.field001[0,*]) flask = reform(data.field001[ii,*]) for t = 0, nco2_ext-1 do begin chk = where(dates ge co2_time_ext[t]-1./24. and dates lt co2_time_ext[t]+1./24.,nchk) if nchk gt 0 then flask_data_ext[t,3,i] = mean(flask[chk]) endfor gd = where(finite(flask_data_ext[*,0,i])) smo_yrsm = lowpass_cfc(flask_data_ext[gd,0,i], BOX=12, EDGE_PFCAST=1) deriv_yr = 12*sg_smooth(smo_yrsm, NLEFT=12/2, NRIGHT=12/2-1, DERIV=1, DELTA=1.0, EDGE_PFCAST=1) smo_gr = deriv_yr / smo_yrsm smo_gr = mean(smo_gr[0:23]) for y = 0, 3 do begin gdy = where(finite(flask_data_ext[*,y,i])) for t = gdy[0]-1, 0, -1 do begin if ci_vars_dc8[i] eq 'CFC-114_WAS' and co2_time_ext[t] le 1998 then smo_gr_adj = 0.02 else smo_gr_adj = smo_gr if ci_vars_dc8[i] eq 'C2Cl4_WAS' then smo_gr_adj *= 3. if ci_vars_dc8[i] eq 'CHCl3_WAS' then smo_gr_adj = 0. flask_data_ext[t,y,i] = (1.-smo_gr_adj/12.) * flask_data_ext[t+1,y,i] endfor endfor flask_data_ext[*,4,i] = flask_data_ext[*,1,i] ; Tropics flask_data_ext[*,5,i] = flask_data_ext[*,2,i] ; NH if ci_vars_dc8[i] eq 'CFC-114_WAS' then flask_data_ext[*,*,i] *= 0.972;3 endif endfor ; -------------- Source latitude distributions. ----------------------------------------------------------------------- ny_source = 65 npeaks_source = 8 nwidths_source = 1 tr_source_lats = findgen(ny_source)-30. nh_source_lats = findgen(ny_source)+12. tr_peaks = findgen(npeaks_source)*4. - 4. nh_peaks = findgen(npeaks_source)*4. + 26. ;peak_widths = [15.,50.] peak_widths = [20.] tr_gauss = fltarr(ny_source,npeaks_source,nwidths_source) & nh_gauss = tr_gauss for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin tr_gauss[*,pp,wi] = exp(-(tr_source_lats-tr_peaks[pp])^2. / (2.*peak_widths[wi])) tr_gauss[*,pp,wi] /= total(tr_gauss[*,pp,wi]) nh_gauss[*,pp,wi] = exp(-(nh_source_lats-nh_peaks[pp])^2. / (2.*peak_widths[wi])) nh_gauss[*,pp,wi] /= total(nh_gauss[*,pp,wi]) endfor ;p = plot(indgen(2),/nodata,xrange=[-15,65],yrange=[0,0.11],xtitle='Latitude',ytitle='Frequency',title='Source Latitude Distributions',dimensions=[700,500],margin=0.1,font_size=11) ;for pp = 0, npeaks_source-1 do p = plot(tr_source_lats,tr_gauss[*,pp,0],thick=2,color='lime green',/overplot) ;for pp = 0, npeaks_source-1 do p = plot(nh_source_lats,nh_gauss[*,pp,0],thick=2,color='medium orchid',/overplot) ;;for wi = 0, nwidths_source-1 do p = plot(tr_source_lats,tr_gauss[*,5,wi],thick=2,color='medium orchid',/overplot) ;;for wi = 0, nwidths_source-1 do p = plot(nh_source_lats,nh_gauss[*,5,wi],thick=2,color='lime green',/overplot) ;txt = text(0.3,0.8,'Tropics',color='lime green',/norm) ;txt = text(0.6,0.8,'N. America',color='medium orchid',/norm) ;p.save,dir+'Plots/SEAC4RS/Source_lat_distributions.png' ; Interpolate surface data to one degree latitude resolution. ny_surf_interp = 117 lat_surf_interp = findgen(ny_surf_interp)-30. surface_data_interp = fltarr(nco2_ext,ny_surf_interp,nciv) for i = 1, nciv-1 do begin if i le 2 then begin slats = [-14.3,19.5,25.7,45.9] yi = interpol(indgen(4),slats,lat_surf_interp) for t = 0, nco2_ext-1 do surface_data_interp[t,*,i] = interpolate(reform(flask_data_ext[t,0:3,i]),yi) endif if (i ge 3 and i le 6) or (i ge 9 and i le 19) then begin slats = [-10.0,10.0,32.5,52.5] yi = interpol(indgen(4),slats,lat_surf_interp) for t = 0, nco2_ext-1 do surface_data_interp[t,*,i] = interpolate(reform(flask_data_ext[t,0:3,i]),yi) endif if i eq 7 or i eq 8 or i eq 20 then begin slats = [-14.3,13.0,41.0,53.0] yi = interpol(indgen(4),slats,lat_surf_interp) for t = 0, nco2_ext-1 do surface_data_interp[t,*,i] = interpolate(reform(flask_data_ext[t,0:3,i]),yi) endif endfor ; Find the tropical and NH surface time series for all of the source region combinations. y1 = where(lat_surf_interp eq nh_source_lats[0]) surf = fltarr(nco2_ext,npeaks_source,nwidths_source,2,nciv) for i = 0, nciv-1 do for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do for t = 0, nco2_ext-1 do begin surf[t,pp,wi,0,i] = total(tr_gauss[*,pp,wi] * reform(surface_data_interp[t,0:ny_source-1,i])) surf[t,pp,wi,1,i] = total(nh_gauss[*,pp,wi] * reform(surface_data_interp[t,y1:y1+ny_source-1,i])) endfor co2_surface_interp = fltarr(nco2_ext,ny_surf_interp) yi = interpol(indgen(ny_ct),lat_ct,lat_surf_interp) for t = 0, nco2_ext-1 do co2_surface_interp[t,*] = interpolate(reform(co2_ct_combo_zm[t,*]),yi) co2_surf = fltarr(nco2_ext,npeaks_source,nwidths_source,2) for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do for t = 0, nco2_ext-1 do begin co2_surf[t,pp,wi,0] = total(tr_gauss[*,pp,wi] * reform(co2_surface_interp[t,0:ny_source-1])) co2_surf[t,pp,wi,1] = total(nh_gauss[*,pp,wi] * reform(co2_surface_interp[t,y1:y1+ny_source-1])) endfor ; Age spectrum calculation. -------------------------------------------------------------------------------------------- ni = 15 ; Number of iterations for minimization loop. ci_lifetimes_theta_profile = fltarr(ni,nciv,nzth) & ci_lifetimes_linear_theta_profile = ci_lifetimes_theta_profile tau_profile = fltarr(nciv,nzth) & tau_profile_source = tau_profile for i = 0, nciv-1 do begin ci_lifetimes_theta_profile[*,i,*] = replicate(ci_lifetimes[i],ni,nzth) ci_lifetimes_linear_theta_profile[*,i,*] = replicate(ci_lifetimes[i],ni,nzth) tau_profile[i,*] = replicate(ci_lifetimes[i],nzth) tau_profile_source[i,*] = replicate(ci_lifetimes[i],nzth) endfor ngt = 121. ngtl = 8.4e3 ngt_1yr = 3.75e3 ngt_tr = 5.1e3 H = 7600. z = 1.7e4 nr = 20 ng = 140 grn_mean_ages = findgen(ng)*2+1. ratios = findgen(nr)*0.1+0.1 Kdiff = [10^(findgen(ng-10)/(ng/3.))+1.,findgen(5)*5e2+1e3,findgen(5)*1e3+4e3] Kdiff_1yr = 10^(findgen(ng)/(ng/3.)+0.5)+100. Kdiff_tr = 10^(findgen(ng)/(ng/2.5))+1. ci_theta_grid = replicate(!values.f_nan,2,nzth,nciv) & cc = replicate(!values.f_nan,nzth,nciv) g_linear_theta_all = fltarr(ngtl,ng,nzth) & gr_linear_theta = dblarr(ngtl,ni,nzth) & trop_frac_g_theta = dblarr(ngt,ng,nzth) trop_frac_gr_linear_theta = gr_linear_theta & trop_fracs_linear_theta_fit = replicate(!values.f_nan,ngtl,ni,nzth) trop_frac_g_linear_theta = dblarr(ngtl,ng,nzth) trop_frac_gr_theta_nh = dblarr(ngt_1yr,ng,nzth) trop_frac_gr_theta_tr = dblarr(ngt_tr,ng,nzth) & bl_frac_gr_theta_region = replicate(!values.f_nan,ng,4,nzth,nciv) & df = fltarr(ng,nzth) & tfk = fltarr(ngtl,ng,nzth) tf_gradj = replicate(!values.f_nan,ng,nzth,nciv) & tfk_fit = tfk & g_theta_best = fltarr(ngtl,nzth) & tfk_best = g_theta_best & tf_gradj_best = cc & bl_gradj = tf_gradj tau_k = tf_gradj & co2_g = fltarr(ny_ct,ng,nzth) & diff_min_i = fltarr(nzth) & sspb = df & seas = fltarr(ngtl,12) co2_g_best_seas = replicate(!values.f_nan,nzth,12) & bl_surf = replicate(!values.f_nan,npeaks_source,nwidths_source,npeaks_source,nwidths_source,ng,nzth,nciv) & tf_surf = bl_surf diffs_source = tf_surf & diffs_source_tot = replicate(!values.f_nan,npeaks_source,nwidths_source,npeaks_source,nwidths_source,ng,nzth) & sspb_source = diffs_source_tot si_co2_ki = replicate(!values.f_nan,npeaks_source^2*nwidths_source^2,ng,nzth) & si_diffs_ki = si_co2_ki & si_sspb_ki = si_co2_ki & sspb_ki_med = df & nh_width = diff_min_i zz = replicate(!values.f_nan,nzth) & tr_width = zz & nh_peak = zz & tr_peak = zz & ki_prof = zz & tfk_best_source = g_theta_best & tf_best_source = cc & tr_best_source_frac = diff_min_i g_best_source = g_theta_best & mean_age_source = zz & modal_age_source = zz & g_best_source_tr = g_best_source & g_best_source_nh = g_best_source mean_age_source_tr = zz & mean_age_source_nh = zz & modal_age_source_tr = zz & modal_age_source_nh = zz & bl_best_source = cc & co2_best_source = diff_min_i mean_age_g_linear_theta_all = df & modal_age_g_linear_theta_all = df & diff_tot_source = diffs_source_tot & bl_nh_best_source = cc & bl_tr_best_source = cc life_grid_linear = [findgen(8e3)*0.5+0.5,findgen(399)*182.5+11.*365.,1e5] life_grid_linear_log = alog10(life_grid_linear) alt_vname = 'GPS_ALT_NASDAT' for z = 0, nzth-1 do begin z_theta = (1e3*reform(er2_was_theta_grid[1,z,where(varnames_er2_was eq alt_vname)]))[0] for ki = 0, ng-1 do begin g_linear_theta_all[*,ki,z] = z_theta/(2.*sqrt(!pi*Kdiff[ki]*life_grid_linear^3)) * exp(z_theta/(2.*H) - Kdiff[ki]*life_grid_linear/(4.*H^2) - z_theta/(4.*Kdiff[ki]*life_grid_linear)) g_linear_theta_all[*,ki,z] /= total(g_linear_theta_all[*,ki,z]) mean_age_g_linear_theta_all[ki,z] = total(life_grid_linear*g_linear_theta_all[*,ki,z]) mm = max(g_linear_theta_all[*,ki,z],mi) modal_age_g_linear_theta_all[ki,z] = life_grid_linear[mi] endfor endfor ;print,mean_age_g_linear_theta_all[*,6],modal_age_g_linear_theta_all[*,6] ; BL fractions for predefined spectra shapes. ;for z = 5, 18 do for ki = 0, ng-1 do for i = 0, ngtl-1 do tfk[i,ki,z] = total(exp(-life_grid_linear/life_grid_linear[i]) * g_linear_theta_all[*,ki,z]) ;save,life_grid_linear,g_linear_theta_all,tfk,filename=dir+'Aircraft/Missions/seac4rs/age_spectra_predef_tfk.sav' restore,dir+'Aircraft/Missions/seac4rs/age_spectra_predef_tfk.sav' zti = where(er2_was['THETA'] ge 350,nm) er2_yrfrac = 2013. + er2_was['Fractional_Day',zti]/365. er2_yrfrac_mean = mean(er2_yrfrac) ;print,er2_yrfrac_mean ti = interpol(indgen(nco2_ext),co2_time_ext,er2_yrfrac_mean-life_grid_linear/365.) tii = where(life_grid_linear le 60) tii2 = where(life_grid_linear gt 60 and life_grid_linear lt 1.09e4) ti_30yr = where(life_grid_linear lt 1.09e4,ngt_30yr) ;vname = 'HFC-134a_WAS' ;i = where(varnames_er2_was eq vname) ;i2 = where(varnames_dc8_was eq vname) ;i3 = where(ci_vars_dc8 eq vname) ;p = plot(indgen(2),/nodata,title=vname+' 340-360K') ;p = plot(er2_was_theta_grid_hist_arr[*,5,i],er2_was_theta_grid_hist[*,5,i],color='blue',thick=3,/overplot) ;p = plot(replicate(er2_was_theta_grid[1,6,i],2),[0.1,0.1],symbol='o',/sym_filled,color='blue',/overplot) ;p = plot(dc8_was_theta_grid_hist_arr[*,5,i2],dc8_was_theta_grid_hist[*,5,i2],color='sky blue',thick=3,/overplot) ;p = plot(dc8_was_theta_grid[1,5:6,i2],[0.1,0.1],symbol='o',/sym_filled,color='sky blue',/overplot) ;p = plot(flask_data_ext[ti[0],1:5,i3],replicate(0.05,5),symbol='o',/sym_filled,color='red',/overplot) ;p = plot(0.97*flask_data_ext[ti[0],1:5,i3],replicate(0.05,5),symbol='o',color='red',/overplot) ;txt = text(0.2,0.8,'ER-2',color='blue',/norm) ;txt = text(0.2,0.76,'DC-8',color='sky blue',/norm) ;print,co2_time_ext[ti[0]],reform(flask_data_ext[ti[0],1:5,i3]) ; Time weighting of age spectra from tropics and NH. day1 = 1. day2 = [30.,40.,50.,60.,80.,100.,120.,150.,200.,250.] ;p = plot(indgen(2),/nodata,xrange=[0,120],xtitle='Age (days)',ytitle='Fraction',title='Partitioning of NH vs. Tropical Age Spectra as a Function of Age',dimensions=[700,500],margin=0.09,font_size=11) ;for di = 0, n_elements(day2)-1 do begin ; t1 = where(life_grid_linear le day1,nt1) ; t2 = where(life_grid_linear gt day1 and life_grid_linear lt day2[di],nt2) ; t3 = where(life_grid_linear ge day2[di],nt3) ; linear1 = 1. - findgen(nt2)/nt2 ; gauss1 = exp(-(life_grid_linear[t2]-day1)^2. / (day2[di]/30.*(day2[di]-day1)*10.)) ; gauss1 -= gauss1[-1] ; ss = gauss1 ; weight_tr = [fltarr(nt1),reverse(ss),fltarr(nt3)+1] ; weight_nh = [fltarr(nt1)+1,ss,fltarr(nt3)] ; tot = weight_tr + weight_nh ; weight_tr /= tot ; weight_nh /= tot ; ; p = plot(life_grid_linear,weight_tr,color='lime green',thick=2,linestyle=(di mod 6),/overplot) ; p = plot(life_grid_linear,weight_nh,color='medium orchid',thick=2,linestyle=(di mod 6),/overplot) ;endfor shp_i = 2 t1 = where(life_grid_linear le day1,nt1) t2 = where(life_grid_linear gt day1 and life_grid_linear lt day2[shp_i],nt2) t3 = where(life_grid_linear ge day2[shp_i],nt3) linear1 = 1. - findgen(nt2)/nt2 gauss1 = exp(-(life_grid_linear[t2]-day1)^2. / (day2[shp_i]/30.*(day2[shp_i]-day1)*10.)) gauss1 -= gauss1[-1] ss = gauss1 weight_tr = [fltarr(nt1),reverse(ss),fltarr(nt3)+1] weight_nh = [fltarr(nt1)+1,ss,fltarr(nt3)] tot = weight_tr + weight_nh weight_tr /= tot weight_nh /= tot shp_i_opt = intarr(nzth) shp_i_opt[6:8] = replicate(2,3) shp_i_opt[9] = 3 shp_i_opt[10] = 5 shp_i_opt[11:17] = replicate(7,7) shp_i_opt[18] = 8 ;p = plot(indgen(2),/nodata,xrange=[0,60],xtickvalues=[0,25,50],xtickname=['0','$\itt_{f}$ /2','$\itt_{f}$'],xminor=0,xtitle='Age',ytitle='Fraction of Total Age Spectra', $ ; title='Partitioning of NH vs. Tropical Age Spectra as a Function of Age',dimensions=[600,450],margin=0.11,font_size=12) ;p = plot([25,25],[0,1],linestyle=1,/overplot) ;p = plot([50,50],[0,1],linestyle=1,/overplot) ;tik = where(life_grid_linear le 60) ;poly = polygon([0,60,reverse(life_grid_linear[tik])],[0,0,reverse(weight_tr[tik])],/data,/fill_background,fill_color='lime green',fill_transparency=80,color='lime green',thick=2) ;poly = polygon([0,life_grid_linear[tik],60,0],[0,weight_tr[tik],1,1],/data,/fill_background,fill_color='medium orchid',fill_transparency=80,color='medium orchid',thick=2) ;p = plot(life_grid_linear,weight_tr,thick=2,/overplot) ;t = text(0.18,0.7,'From N. American',font_size=11,color='medium orchid',/norm) ;t = text(0.18,0.66,'Surface',font_size=11,color='medium orchid',/norm) ;t = text(0.52,0.3,'From Tropical',font_size=11,color='lime green',/norm) ;t = text(0.52,0.26,'Surface',font_size=11,color='lime green',/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_spectra_partitioning.png' ;p = plot(indgen(2),/nodata,xrange=[0,120]) ;p = plot(life_grid_linear,weight_tr,color='blue',thick=2,/overplot) ;p = plot(life_grid_linear,weight_nh,color='orange',thick=2,/overplot) ; Weighting shape. shapes = ['gauss_30day','gauss_40day','gauss_50day','gauss_60day','gauss_80day','gauss_100day','gauss_120day','gauss_150day','gauss_200day','gauss_250day'] print,shapes[shp_i] ; BL CO2 for predefined spectra shapes. ;co2_g_all = fltarr(npeaks_source,nwidths_source,npeaks_source,nwidths_source,ng,nzth) & co2_g_all_diffs = co2_g_all ;for z = 5, 18 do for ki = 0, ng-1 do begin ; gr_theta_30yr = g_linear_theta_all[ti_30yr,ki,z] ; gr_theta_30yr /= total(gr_theta_30yr) ; gr_sub1 = weight_nh * gr_theta_30yr ; gr_sub2 = weight_tr * gr_theta_30yr ; for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin ; co2_tr_int = interpolate(co2_surf[*,pp,wi,0],ti) ; for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do begin ; co2_nh_int = interpolate(co2_surf[*,pp2,wi2,1],ti) ; co2_g_all[pp,wi,pp2,wi2,ki,z] = total(gr_sub1 * co2_nh_int) + total(gr_sub2 * co2_tr_int) ; endfor ; endfor ; co2_g_all_diffs[*,*,*,*,ki,z] = co2_g_all[*,*,*,*,ki,z] - (er2_was_theta_grid[1,z,where(varnames_er2_was eq 'CO2_HUPCRS')])[0] ; if ki eq 0 then print,z ;endfor ;save,weight_nh,weight_tr,co2_g_all,co2_g_all_diffs,filename=dir+'Aircraft/Missions/seac4rs/co2_predef_'+shapes[shp_i]+'.sav' restore,dir+'Aircraft/Missions/seac4rs/co2_predef_'+shapes[shp_i]+'.sav' ; BL tracer values for predefined spectra shapes. ;bl_g_all = fltarr(npeaks_source,nwidths_source,npeaks_source,nwidths_source,ng,nzth,nciv) & bl_g_nh = bl_g_all & bl_g_tr = bl_g_all ;for i = 1, nciv-1 do for z = 5, 18 do for ki = 0, ng-1 do begin ; gr_theta_30yr = g_linear_theta_all[ti_30yr,ki,z] ; gr_theta_30yr /= total(gr_theta_30yr) ; gr_sub1 = weight_nh * gr_theta_30yr ; gr_sub2 = weight_tr * gr_theta_30yr ; for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin ; bl_tr_int = interpolate(surf[*,pp,wi,0,i],ti) ; for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do begin ; bl_nh_int = interpolate(surf[*,pp2,wi2,1,i],ti) ; bl_g_nh[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub1 * bl_nh_int) ; bl_g_tr[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub2 * bl_tr_int) ; bl_g_all[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub1 * bl_nh_int) + total(gr_sub2 * bl_tr_int) ; endfor ; endfor ; if z eq 6 and ki eq 0 then print,i ;endfor ;save,bl_g_all,bl_g_nh,bl_g_tr,filename=dir+'Aircraft/Missions/seac4rs/bl_predef_'+shapes[shp_i]+'.sav' restore,dir+'Aircraft/Missions/seac4rs/bl_predef_'+shapes[shp_i]+'.sav' for i = 0, nciv-1 do begin er2i = where(strmatch(strmid(varnames_er2_was,0,strlen(ci_vars_er2[i])),ci_vars_er2[i],/fold_case) eq 1) dc8i = where(strmatch(strmid(varnames_dc8_was,0,strlen(ci_vars_dc8[i])),ci_vars_dc8[i],/fold_case) eq 1) ; Boundary layer values. case ci_vars_dc8[i] of 'n-Butane_WAS' : begin er2_was_theta_grid[*,11,er2i] = !values.f_nan end else : endcase if i eq 3 or i eq 6 then ii = [0,2] else ii = [1,2] ci_theta_grid[*,0:4,i] = dc8_was_theta_grid[ii,0:4,dc8i] ci_theta_grid[*,5:-1,i] = er2_was_theta_grid[ii,5:-1,er2i] if i eq 9 then ci_theta_grid[0,6:-1,i] = smooth(ci_theta_grid[0,6:-1,i],3,/edge_truncate) if i eq 13 then for gg = 0, 2 do ci_theta_grid[0,6:-3,i] = smooth(ci_theta_grid[0,6:-3,i],3,/edge_truncate) if i eq 15 then for gg = 0, 5 do ci_theta_grid[0,6:-3,i] = smooth(ci_theta_grid[0,6:-3,i],3,/edge_truncate) endfor ; Restore the CLaMS spectra. restore,dir+'Models/CLaMS/CLaMS_age_spec_SEAC4RS.sav' ti_c = interpol(indgen(nst_ext),365.*yr_spec_ext,life_grid_linear) zi_c = interpol(indgen(nth_spec),theta_spec,theta_grid) age_spec_tr_clams_seac4rs_interp = interpolate(age_spec_tr_clams_seac4rs,ti_c,zi_c,/grid) age_spec_nh_clams_seac4rs_interp = interpolate(age_spec_nh_clams_seac4rs,ti_c,zi_c,/grid) age_spec_tot_clams_seac4rs_interp = age_spec_tr_clams_seac4rs_interp tf_clams = fltarr(ngtl,nzth) for z = 6, nzth-1 do begin age_spec_tot_clams_seac4rs_interp[*,z] = age_spec_tr_clams_seac4rs_interp[*,z] + age_spec_nh_clams_seac4rs_interp[*,z] frac_tr = age_spec_tr_clams_seac4rs_interp[*,z] / age_spec_tot_clams_seac4rs_interp[*,z] chk = where(~finite(frac_tr),nchk) if nchk gt 0 then frac_tr[chk] = 0. age_spec_tot_clams_seac4rs_interp[*,z] /= total(age_spec_tot_clams_seac4rs_interp[*,z]) age_spec_tr_clams_seac4rs_interp[*,z] = frac_tr * age_spec_tot_clams_seac4rs_interp[*,z] age_spec_nh_clams_seac4rs_interp[*,z] = (1.-frac_tr) * age_spec_tot_clams_seac4rs_interp[*,z] for i = 0, ngtl-1 do tf_clams[i,z] = total(exp(-life_grid_linear/life_grid_linear[i]) * age_spec_tot_clams_seac4rs_interp[*,z]) endfor bl_clams_all = fltarr(npeaks_source,nwidths_source,npeaks_source,nwidths_source,nzth,nciv) & tf_clams_all = bl_clams_all co2_clams_all = fltarr(npeaks_source,nwidths_source,npeaks_source,nwidths_source,nzth) for z = 6, 18 do begin age_spec_tr_30yr = age_spec_tr_clams_seac4rs_interp[ti_30yr,z] age_spec_nh_30yr = age_spec_nh_clams_seac4rs_interp[ti_30yr,z] age_spec_tot_30yr = age_spec_tot_clams_seac4rs_interp[ti_30yr,z] age_spec_tr_30yr /= total(age_spec_tot_30yr) age_spec_nh_30yr /= total(age_spec_tot_30yr) ; if z eq 7 then print,(total(age_spec_tr_clams_seac4rs_interp[*,z]) + total(age_spec_nh_clams_seac4rs_interp[*,z])),(total(age_spec_tr_30yr) + total(age_spec_nh_30yr)), $ ; total(age_spec_tot_clams_seac4rs_interp[*,z]) for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do begin for i = 1, nciv-1 do begin bl_tr_int = interpolate(surf[*,pp,wi,0,i],ti) bl_nh_int = interpolate(surf[*,pp2,wi2,1,i],ti) bl_clams_all[pp,wi,pp2,wi2,z,i] = total(age_spec_nh_30yr * bl_nh_int) + total(age_spec_tr_30yr * bl_tr_int) endfor co2_tr_int = interpolate(co2_surf[*,pp,wi,0],ti) co2_nh_int = interpolate(co2_surf[*,pp2,wi2,1],ti) co2_clams_all[pp,wi,pp2,wi2,z] = total(age_spec_nh_30yr * co2_nh_int) + total(age_spec_tr_30yr * co2_tr_int) endfor for i = 1, nciv-1 do tf_clams_all[*,*,*,*,z,i] = ci_theta_grid[0,z,i]/bl_clams_all[*,*,*,*,z,i] endfor chk = where(tf_clams_all gt 1,nchk) if nchk gt 1 then tf_clams_all[chk] = 1.0 ; ------------------------------------------------------------------------------------------------------------------------------------ ; Begin calculation with theta average profiles. ; ------------------------------------------------------------------------------------------------------------------------------------ restore,filename=dir+'Aircraft/Missions/seac4rs/theta_avg_profiles.sav' ci_theta_grid[0,*,*] = ci_theta_grid_ind_adj ; -------------------------------------- Age spectra for theta grid source regions. --------------------------------------------------------------------------------------------------- doit = 1 ; *********************************** if doit then begin for z = 6, 18 do begin restore,dir+'Aircraft/Missions/seac4rs/co2_predef_'+shapes[shp_i_opt[z]]+'.sav' restore,dir+'Aircraft/Missions/seac4rs/bl_predef_'+shapes[shp_i_opt[z]]+'.sav' gd = where(finite(reform(ci_theta_grid[0,z,1:-1])),ngd) if z ge 7 then begin tau_profile_source[*,z] = replicate(!values.f_nan,nciv) tau_profile_source[gd+1,z] = tau_profile_source[gd+1,z-1] endif tau = tau_profile_source[gd+1,z] tau_i = interpol(indgen(ngtl),life_grid_linear,tau) for ki = 0, ng-1 do begin for i = 1, nciv-1 do begin tf_surf[*,*,*,*,ki,z,i] = ci_theta_grid[0,z,i]/bl_g_all[*,*,*,*,ki,z,i] temp = tf_surf[*,*,*,*,ki,z,i] chk = where(temp gt 1,nchk) if nchk gt 1 then temp[chk] = 1.0 tf_surf[*,*,*,*,ki,z,i] = temp endfor ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,z],tau_i) for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin diffs_source[pp,wi,pp2,wi2,ki,z,gd+1] = (tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] - tfk_i)^2 diffs_source_tot[pp,wi,pp2,wi2,ki,z] = total(diffs_source[pp,wi,pp2,wi2,ki,z,gd+1]) median_ratio = median(alog(tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] / tfk_i)) sspb_source[pp,wi,pp2,wi2,ki,z] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor diff_tot_source[*,*,*,*,ki,z] = 10. * abs(co2_g_all_diffs[*,*,*,*,ki,z]) + abs(sspb_source[*,*,*,*,ki,z]) si_co2_ki[*,ki,z] = sort(abs(co2_g_all_diffs[*,*,*,*,ki,z])) si_diffs_ki[*,ki,z] = sort(diffs_source_tot[*,*,*,*,ki,z]) si_sspb_ki[*,ki,z] = sort(sspb_source[*,*,*,*,ki,z]) sspb_ki_med[ki,z] = median(sspb_source[*,*,*,*,ki,z]) co2_temp = co2_g_all_diffs[*,*,*,*,ki,z] diffs_temp = diffs_source_tot[*,*,*,*,ki,z] sspb_temp = sspb_source[*,*,*,*,ki,z] ; if z eq 6 then print,co2_temp[si_co2_ki[0:5,ki,z]],sspb_temp[si_sspb_ki[0:5,ki,z]],sspb_ki_med[ki,z] endfor ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2*nwidths_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min if z le 10 then kd = 3 else kd = 2 if z eq 6 then ki_max = ng-1 else ki_max = ki_prof[z-1] - kd for ki = 0, ki_max do begin diff_tot = diff_tot_source[*,*,*,*,ki,z] nh_width_tot = fix(indgen(nsources)/(npeaks_source^2*nwidths_source)) nh_peak_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) tr_width_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source - nh_peak_tot*npeaks_source*nwidths_source) / (npeaks_source)) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) si0 = 0 if z ge 7 then begin si = sort(diff_tot) trp = fix(tr_peak[z-1]) si1 = si[0] chk = where(tr_peak_tot ge trp-1 and tr_peak_tot le trp) d_tot_si = sort(diff_tot[chk]) si0 = d_tot_si[0] ; if z eq 8 and ki ge 106 and ki le 110 then print,ki,diff_tot[si[0:10]],diff_tot[chk[d_tot_si[0:10]]],si0,si1,tr_peak_tot[chk[d_tot_si[0:10]]] diff_val[ki] = diff_tot[chk[si0]] nh_w[ki] = nh_width_tot[chk[si0]] nh_p[ki] = nh_peak_tot[chk[si0]] tr_w[ki] = tr_width_tot[chk[si0]] tr_p[ki] = tr_peak_tot[chk[si0]] endif else begin d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[si0]] nh_w[ki] = nh_width_tot[d_tot_si[si0]] nh_p[ki] = nh_peak_tot[d_tot_si[si0]] tr_w[ki] = tr_width_tot[d_tot_si[si0]] tr_p[ki] = tr_peak_tot[d_tot_si[si0]] endelse diff_min[ki] = d_tot_si[si0] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_prof[z] = ki_min ; co2_d = abs(co2_g_all_diffs[*,*,*,*,ki_min,z]) ; co2_si = si_co2_ki[*,ki_min,z] ; co2_diffs = abs(co2_g_all_diffs[*,*,*,*,ki_min,z]) ; gd = where(co2_diffs[co2_si] le 0.2) ; co2_i = co2_si[gd] ; diff_tot = 10.*co2_d + abs(sspb) ; d_tot_si = sort(diff_tot) ; nh_width_tot = fix(d_tot_si[0:10]/(npeaks_source^2*nwidths_source)) ; nh_peak_tot = fix((d_tot_si[0:10] - nh_width_tot*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) ; tr_width_tot = fix((d_tot_si[0:10] - nh_width_tot*npeaks_source^2*nwidths_source - nh_peak_tot*npeaks_source*nwidths_source) / (npeaks_source)) ; tr_peak_tot = fix(d_tot_si[0:10] mod npeaks_source) ; if z eq 8 then print,ki_min,diff_tot[d_tot_si[0:10]],tr_peak_tot,tr_width_tot,nh_peak_tot,nh_width_tot ; nh_width_co2 = fix(co2_i/(npeaks_source^2*nwidths_source)) ; nh_peak_co2 = fix((co2_i - nh_width_co2*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) ; tr_width_co2 = fix((co2_i - nh_width_co2*npeaks_source^2*nwidths_source - nh_peak_co2*npeaks_source*nwidths_source) / (npeaks_source)) ; tr_peak_co2 = fix(co2_i mod npeaks_source) ; if z eq 6 then print,co2_diffs[co2_si[gd]],tr_peak_co2,tr_width_co2,nh_peak_co2,nh_width_co2 ; nh_width_sspb = fix(sspb_i/(npeaks_source^2*nwidths_source)) ; nh_peak_sspb = fix((sspb_i - nh_width_sspb*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) ; tr_width_sspb = fix((sspb_i - nh_width_sspb*npeaks_source^2*nwidths_source - nh_peak_sspb*npeaks_source*nwidths_source) / (npeaks_source)) ; tr_peak_sspb = fix(sspb_i mod npeaks_source) ; if z eq 6 then print,tr_peak_sspb,tr_width_sspb,nh_peak_sspb,nh_width_sspb tr_peak[z] = tr_p[ki_min] tr_width[z] = tr_w[ki_min] nh_peak[z] = nh_p[ki_min] nh_width[z] = nh_w[ki_min] g_best_source[*,z] = g_linear_theta_all[*,ki_min,z] g_best_source_tr[*,z] = weight_tr * g_best_source[*,z] g_best_source_nh[*,z] = weight_nh * g_best_source[*,z] tr_best_source_frac[z] = total(g_best_source_tr[*,z]) / total(g_best_source[*,z]) ; Adjust the lifetimes to best fit the spectra. bl_best_source[z,*] = bl_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_nh_best_source[z,*] = bl_g_nh[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_tr_best_source[z,*] = bl_g_tr[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] tfk_best_source[*,z] = tfk[*,ki_min,z] tf_best_source[z,*] = tf_surf[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] co2_best_source[z] = co2_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z] tfi = interpol(indgen(ngtl),tfk_best_source[*,z],tf_best_source[z,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_profile_source[i,z] = tau endif endfor ; if z eq 6 then print,tau_profile_source[*,z] mean_age_source[z] = total(life_grid_linear*g_best_source[*,z]) mean_age_source_tr[z] = total(life_grid_linear*g_best_source_tr[*,z]) mean_age_source_nh[z] = total(life_grid_linear*g_best_source_nh[*,z]) mm = max(g_best_source[*,z],mi) modal_age_source[z] = life_grid_linear[mi] mm = max(g_best_source_tr[*,z],mi) modal_age_source_tr[z] = life_grid_linear[mi] mm = max(g_best_source_nh[*,z],mi) modal_age_source_nh[z] = life_grid_linear[mi] endfor tau_profile_source_init = tau_profile_source print,ki_prof[5:18],tr_peak[5:18],nh_peak[5:18] endif ; ---------------------------------------------------------------------------------------------------------------------------------- ; Compute age spectra for each individual WAS measurement above 340K, primarily to identify pollution source outliers in profiles. ; ---------------------------------------------------------------------------------------------------------------------------------- theta_ind = er2_was['THETA',zti] zi = interpol(indgen(nzth),theta_grid,theta_ind) alt_ind = er2_was[alt_vname,zti] lon_ind = er2_was['LONGITUDE',zti] lat_ind = er2_was['LATITUDE',zti] tpause_ind = er2_was['TropopausePALT1_MTP',zti] tpause_theta_ind = er2_was['TropopausePotTemp1_MTP',zti] co2_ind = er2_was['CO2_HUPCRS',zti] o3_ind = er2_was['O3_UAS',zti] chk = where(co2_ind gt 0,nchk) print,nm,nchk ;bl_g_ind = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm,nciv) & bl_g_nh_ind = bl_g_ind ;for aa = 0, nm-1 do begin ; for i = 1, nciv-1 do for ki = 0, ng-1 do for pp = 0, npeaks_source-1 do for pp2 = 0, npeaks_source-1 do begin ; bl_g_ind[pp,pp2,ki,aa,i] = interpolate(bl_g_all[pp,0,pp2,0,ki,*,i],zi[aa]) ; bl_g_nh_ind[pp,pp2,ki,aa,i] = interpolate(bl_g_nh[pp,0,pp2,0,ki,*,i],zi[aa]) ; endfor ;endfor ;save,bl_g_ind,bl_g_nh_ind,filename=dir+'Aircraft/Missions/seac4rs/bl_predef_ind_'+shapes[shp_i]+'.sav' nh = 5e2 dn = 1e-2 hmax = 10. H_i = findgen(nh)*hmax/nh doit = 1 ; *********************************** Use source region spectra, identify outliers from pollution sources, make initial adjustments to theta average profiles. ************************ if doit then begin mean_prof_rel = replicate(!values.f_nan,nm,nciv) & scale_nh = replicate(1.0,nm,nciv) tf_ind_init = mean_prof_rel & tf_ind_init_mean_rel = mean_prof_rel & hist_mean_rel = fltarr(nh,nzth,nciv) & hist_all_mean_rel = fltarr(nh,nciv) & hist_tf_ind_init_mean_rel = hist_all_mean_rel bl_nh_ind = mean_prof_rel & bl_all_ind = mean_prof_rel & tf_ind_adj = mean_prof_rel & tf_ind_adj_mean_rel = mean_prof_rel & hist_tf_ind_adj_mean_rel = hist_all_mean_rel hist_tf_ind_init_mean_rel_prof = hist_mean_rel & hist_tf_ind_init_prof = hist_mean_rel & tf_ind_adj_mean = replicate(!values.f_nan,nzth,nciv) ci_theta_grid_adj = reform(ci_theta_grid[0,*,*]) & tf_ind_init_mean_rel_tot_low = replicate(!values.f_nan,nm) surface_data = flask_data_ext restore,dir+'Aircraft/Missions/seac4rs/bl_predef_ind_'+shapes[shp_i]+'.sav' for aa = 0, nm-1 do begin zii = zi[aa] if zii gt nzth-1 then zii = nzth-1 if zii lt 6 then zii = 6 for i = 1, nciv-1 do begin if er2_was[ci_vars_er2[i],zti[aa]] gt 0 then begin tf_ind_init[aa,i] = er2_was[ci_vars_er2[i],zti[aa]] / bl_best_source[round(zii),i] if finite(tf_best_source[zii,i]) then tf_ind_init_mean_rel[aa,i] = tf_ind_init[aa,i] / tf_best_source[zii,i] else begin chk = where(finite(tf_best_source[*,i])) mm = min(abs(zii-chk),mi) ; if i eq 4 and aa eq 2 then print,aa,zii,tf_best_source[zi[aa],i],chk[mi],tf_best_source[chk[mi],i],er2_was[ci_vars_er2[i],zti[aa]],tf_ind_init[aa,i],tf_ind_init_mean_rel[aa,i] tf_ind_init_mean_rel[aa,i] = tf_ind_init[aa,i] / tf_best_source[chk[mi],i] endelse ; if i eq 3 and aa eq 8 then print,aa,zii,tf_best_source[zi[aa],i],er2_was[ci_vars_er2[i],zti[aa]],tf_ind_init[aa,i],tf_ind_init_mean_rel[aa,i] chk = where(finite(reform(ci_theta_grid[0,*,i]))) if chk[-1] ge zi[aa] then mean_prof_rel[aa,i] = er2_was[ci_vars_er2[i],zti[aa]] / interpolate(reform(ci_theta_grid[0,*,i]),zii) $ else mean_prof_rel[aa,i] = 1.0 if theta_ind[aa] lt 420 then begin if i le 8 and mean_prof_rel[aa,i] gt 1.3 then surface_data[*,5,i] *= mean_prof_rel[aa,i] if i gt 8 and i le 13 and mean_prof_rel[aa,i] gt 1.03 then surface_data[*,5,i] *= mean_prof_rel[aa,i] if i gt 13 and mean_prof_rel[aa,i] gt 1.01 then surface_data[*,5,i] *= mean_prof_rel[aa,i] endif endif endfor gd = where(finite(reform(tf_ind_init_mean_rel[aa,1:10])),ngd) chk = where(reform(tf_ind_init_mean_rel[aa,gd+1]) lt 1,nchk) if nchk eq ngd then tf_ind_init_mean_rel_tot_low[aa] = 1 else tf_ind_init_mean_rel_tot_low[aa] = 0 ; if aa eq 516 then print,reform(tf_ind_init_mean_rel[aa,1:10]),tf_ind_init_mean_rel_tot_low[aa] endfor tf_ind_adj = tf_ind_init for z = 6, 18 do begin zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) if nn gt 0 then begin for i = 1, nciv-1 do begin H = histogram(mean_prof_rel[zi_th,i],nbins=nh,max=10,min=0) hist_mean_rel[*,z,i] = H chk = where(finite(tf_ind_init_mean_rel[zi_th,i]),ngd) if ngd ge 5 then begin H = histogram(tf_ind_init_mean_rel[zi_th,i],nbins=nh,max=10,min=0) hist_tf_ind_init_mean_rel_prof[*,z,i] = H H = histogram(tf_ind_init[zi_th[chk],i],nbins=nh,max=10,min=0) hist_tf_ind_init_prof[*,z,i] = H tf_ind_adj_mean[z,i] = mean(tf_ind_adj[zi_th,i],/nan) endif endfor endif endfor for i = 1, nciv-1 do begin gd = where(mean_prof_rel[*,i] ne 1) H = histogram(mean_prof_rel[gd,i],nbins=nh,max=10,min=0) hist_all_mean_rel[*,i] = H H = histogram(tf_ind_init_mean_rel[*,i],nbins=nh,max=10,min=0) hist_tf_ind_init_mean_rel[*,i] = H endfor hist_tf_ind_adj_prof = hist_tf_ind_init_prof ; Adjust the pollution source BL fractions. Use elevated propane BL fracs as main indicator. for z = 6, 10 do begin zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) gd = where(finite(tf_ind_adj[zi_th,3]),nng) si = reverse(sort(tf_ind_adj[zi_th[gd],3])) for x = 0, nng-1 do begin for i = 1, nciv-1 do begin chk = where(finite(tf_ind_adj[zi_th,i]),ngd) if ngd gt 1 then begin med = median(tf_ind_adj[zi_th[chk],i]) mm = max(hist_tf_ind_adj_prof[*,z,i],mi) peak_val = H_i[mi] aa = zi_th[gd[si[x]]] zii = zi[aa] if zii gt nzth-1 then zii = nzth-1 if zii lt 6 then zii = 6 mean = interpolate(reform(tf_ind_adj_mean[*,i]),zii) if i le 10 then lim = 1.2 else lim = 1.1 if tf_ind_adj[aa,i] gt lim * mean and tf_ind_init_mean_rel_tot_low[aa] eq 0 then begin ; if (tf_ind_adj[aa,i] gt lim * mean or tf_ind_adj[aa,i] lt 1./lim * mean) and tf_ind_init_mean_rel_tot_low[aa] eq 0 then begin ; Scale up the NH BL values and recalculate the total BL values. scale_nh[aa,i] = tf_ind_adj[aa,i] / med ;mean ;peak_val ; if i ge 3 and i le 8 then scale_nh[aa,i] *= (1./(1.-tr_best_source_frac[zii])) gr_theta_30yr = g_linear_theta_all[ti_30yr,ki_prof[z],z] gr_theta_30yr /= total(gr_theta_30yr) gr_sub1 = weight_nh * gr_theta_30yr bl_tr_ind = bl_tr_best_source[z,i] bl_nh_int = scale_nh[aa,i] * interpolate(surf[*,nh_peak[z],0,1,i],ti) bl_nh_ind[aa,i] = total(gr_sub1 * bl_nh_int) bl_all_ind[aa,i] = bl_tr_ind + bl_nh_ind[aa,i] tf_ind_adj[aa,i] = er2_was[ci_vars_er2[i],zti[aa]] / bl_all_ind[aa,i] chk = where(finite(tf_ind_adj[zi_th,i]),ngd) if ngd ge 5 then begin H = histogram(tf_ind_adj[zi_th[chk],i],nbins=nh,max=10,min=0) hist_tf_ind_adj_prof[*,z,i] = H tf_ind_adj_mean[z,i] = mean(tf_ind_adj[zi_th,i],/nan) endif endif endif endfor endfor for i = 1, nciv-1 do if finite(tf_ind_adj_mean[z,i]) then ci_theta_grid_adj[z,i] *= (tf_ind_adj_mean[z,i] / tf_best_source[z,i]) endfor ; print,tf_ind_adj[8,3],scale_nh[8,3] ; for i = 1, 1 do begin ; ; ; Gaussian fit to BL fracs around the mean profile to identify outliers. ; hsm = smooth(hist_tf_ind_init_mean_rel[*,i],11,/edge_truncate) ; hsm /= max(hsm,mi) ; chk = where(hsm ge 0.5) ; half = 0.1*(H_i[chk[-1]] - H_i[chk[0]]) ; peak_val = H_i[mi] ; gauss = exp(-(H_i-H_i[mi])^2. / (2.*half)) ; gauss /= max(gauss) ; lims = where(gauss ge 0.01) ; hi_vals = H_i[lims[-1]+1:-1] ;; print,peak_val,hi_vals[0] ; ; for aa = 0, nm-1 do begin ; ; zii = zi[aa] ; if zii gt nzth-1 then zii = nzth-1 ; if zii lt 5 then zii = 5 ; ; if finite(tf_ind_init_mean_rel[aa,i]) then begin ; if tf_ind_init_mean_rel[aa,i] ge hi_vals[0] and theta_ind[aa] lt 420 then begin ; ; scale_nh[aa,i] = tf_ind_init_mean_rel[aa,i] / peak_val ; ; ; Scale up the NH BL values and recalculate the total BL values. ; gr_theta_30yr = g_linear_theta_all[ti_30yr,ki_prof[round(zii)],round(zii)] ; gr_theta_30yr /= total(gr_theta_30yr) ; gr_sub1 = weight_nh * gr_theta_30yr ; gr_sub2 = weight_tr * gr_theta_30yr ; bl_tr_ind = bl_tr_best_source[round(zii),i] ; ; bl_nh_int = scale_nh[aa,i] * interpolate(surf[*,nh_peak[round(zii)],0,1,i],ti) ; bl_nh_ind[aa,i] = total(gr_sub1 * bl_nh_int) ; ; bl_nh_ind[aa,i] = scale_nh[aa,i] * bl_nh_best_source[round(zii),i] ; bl_all_ind[aa,i] = bl_tr_ind + bl_nh_ind[aa,i] ; tf_ind_adj[aa,i] = er2_was[ci_vars_er2[i],zti[aa]] / bl_all_ind[aa,i] ; ; print,aa,zi[aa],scale_nh[aa,i],mean_prof_rel[aa,i],bl_nh_best_source[round(zii),i],bl_nh_ind[aa,i],bl_all_ind[aa,i],bl_best_source[round(zii),i],ci_theta_grid[0,round(zii),i], $ ; ; er2_was[ci_vars_er2[i],zti[aa]] ; ; if finite(tf_best_source[zii,i]) then tf_ind_adj_mean_rel[aa,i] = tf_ind_adj[aa,i] / tf_best_source[zii,i] else begin ; chk = where(finite(tf_best_source[*,i])) ; mm = min(abs(zii-chk),mi) ; tf_ind_adj_mean_rel[aa,i] = tf_ind_adj[aa,i] / tf_best_source[chk[mi],i] ; endelse ; ; endif else begin ; tf_ind_adj[aa,i] = tf_ind_init[aa,i] ; if finite(tf_best_source[zii,i]) then tf_ind_adj_mean_rel[aa,i] = tf_ind_adj[aa,i] / tf_best_source[zii,i] else begin ; chk = where(finite(tf_best_source[*,i])) ; mm = min(abs(zii-chk),mi) ; tf_ind_adj_mean_rel[aa,i] = tf_ind_adj[aa,i] / tf_best_source[chk[mi],i] ; endelse ; endelse ; endif ; ; endfor ; ; ; bl_nh_int = interpolate(surf[*,pp2,wi2,1,i],ti) ; ; bl_g_nh[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub1 * bl_nh_int) ; ; bl_g_tr[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub2 * bl_tr_int) ;; bl_g_all[pp,wi,pp2,wi2,ki,z,i] = total(gr_sub1 * bl_nh_int) + total(gr_sub2 * bl_tr_int) ; ; endfor ; ; endfor ; ;endfor ; ; endfor ; ; for i = 1, nciv-1 do begin ; H = histogram(tf_ind_adj_mean_rel[*,i],nbins=nh,max=10,min=0) ; hist_tf_ind_adj_mean_rel[*,i] = H ; endfor endif ; ------------------------------------------------------------------------------------ ; Redo theta average calculation with adjusted mean profiles ; ------------------------------------------------------------------------------------ ; ; -------------------------------------- Age spectra for theta grid source regions. --------------------------------------------------------------------------------------------------- doit = 0 ; *********************************** if doit then begin for z = 6, 18 do begin gd = where(finite(reform(ci_theta_grid_adj[z,1:-1])),ngd) if z ge 7 then begin tau_profile_source[*,z] = replicate(!values.f_nan,nciv) tau_profile_source[gd+1,z] = tau_profile_source[gd+1,z-1] endif tau = tau_profile_source[gd+1,z] tau_i = interpol(indgen(ngtl),life_grid_linear,tau) for ki = 0, ng-1 do begin for i = 1, nciv-1 do begin tf_surf[*,*,*,*,ki,z,i] = ci_theta_grid_adj[z,i]/bl_g_all[*,*,*,*,ki,z,i] temp = tf_surf[*,*,*,*,ki,z,i] chk = where(temp gt 1,nchk) if nchk gt 1 then temp[chk] = 1.0 tf_surf[*,*,*,*,ki,z,i] = temp endfor ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,z],tau_i) for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin diffs_source[pp,wi,pp2,wi2,ki,z,gd+1] = (tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] - tfk_i)^2 diffs_source_tot[pp,wi,pp2,wi2,ki,z] = total(diffs_source[pp,wi,pp2,wi2,ki,z,gd+1]) median_ratio = median(alog(tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] / tfk_i)) sspb_source[pp,wi,pp2,wi2,ki,z] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor diff_tot_source[*,*,*,*,ki,z] = 10. * abs(co2_g_all_diffs[*,*,*,*,ki,z]) + abs(sspb_source[*,*,*,*,ki,z]) si_co2_ki[*,ki,z] = sort(abs(co2_g_all_diffs[*,*,*,*,ki,z])) si_diffs_ki[*,ki,z] = sort(diffs_source_tot[*,*,*,*,ki,z]) si_sspb_ki[*,ki,z] = sort(sspb_source[*,*,*,*,ki,z]) sspb_ki_med[ki,z] = median(sspb_source[*,*,*,*,ki,z]) co2_temp = co2_g_all_diffs[*,*,*,*,ki,z] diffs_temp = diffs_source_tot[*,*,*,*,ki,z] sspb_temp = sspb_source[*,*,*,*,ki,z] ; if z eq 6 then print,co2_temp[si_co2_ki[0:5,ki,z]],sspb_temp[si_sspb_ki[0:5,ki,z]],sspb_ki_med[ki,z] endfor ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2*nwidths_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min if z le 10 then kd = 3 else kd = 2 if z eq 6 then ki_max = ng-1 else ki_max = ki_prof[z-1] - kd for ki = 0, ki_max do begin diff_tot = diff_tot_source[*,*,*,*,ki,z] nh_width_tot = fix(indgen(nsources)/(npeaks_source^2*nwidths_source)) nh_peak_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) tr_width_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source - nh_peak_tot*npeaks_source*nwidths_source) / (npeaks_source)) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) si0 = 0 if z ge 7 then begin si = sort(diff_tot) trp = fix(tr_peak[z-1]) si1 = si[0] chk = where(tr_peak_tot ge trp-1 and tr_peak_tot le trp) d_tot_si = sort(diff_tot[chk]) si0 = d_tot_si[0] diff_val[ki] = diff_tot[chk[si0]] nh_w[ki] = nh_width_tot[chk[si0]] nh_p[ki] = nh_peak_tot[chk[si0]] tr_w[ki] = tr_width_tot[chk[si0]] tr_p[ki] = tr_peak_tot[chk[si0]] endif else begin d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[si0]] nh_w[ki] = nh_width_tot[d_tot_si[si0]] nh_p[ki] = nh_peak_tot[d_tot_si[si0]] tr_w[ki] = tr_width_tot[d_tot_si[si0]] tr_p[ki] = tr_peak_tot[d_tot_si[si0]] endelse diff_min[ki] = d_tot_si[si0] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_prof[z] = ki_min tr_peak[z] = tr_p[ki_min] tr_width[z] = tr_w[ki_min] nh_peak[z] = nh_p[ki_min] nh_width[z] = nh_w[ki_min] g_best_source[*,z] = g_linear_theta_all[*,ki_min,z] g_best_source_tr[*,z] = weight_tr * g_best_source[*,z] g_best_source_nh[*,z] = weight_nh * g_best_source[*,z] ; Adjust the lifetimes to best fit the spectra. bl_best_source[z,*] = bl_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_nh_best_source[z,*] = bl_g_nh[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_tr_best_source[z,*] = bl_g_tr[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] tfk_best_source[*,z] = tfk[*,ki_min,z] tf_best_source[z,*] = tf_surf[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] co2_best_source[z] = co2_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z] tfi = interpol(indgen(ngtl),tfk_best_source[*,z],tf_best_source[z,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_profile_source[i,z] = tau endif endfor ; if z eq 6 then print,tau_profile_source[*,z] mean_age_source[z] = total(life_grid_linear*g_best_source[*,z]) mean_age_source_tr[z] = total(life_grid_linear*g_best_source_tr[*,z]) mean_age_source_nh[z] = total(life_grid_linear*g_best_source_nh[*,z]) mm = max(g_best_source[*,z],mi) modal_age_source[z] = life_grid_linear[mi] mm = max(g_best_source_tr[*,z],mi) modal_age_source_tr[z] = life_grid_linear[mi] mm = max(g_best_source_nh[*,z],mi) modal_age_source_nh[z] = life_grid_linear[mi] endfor print,ki_prof[5:18],tr_peak[5:18],nh_peak[5:18] endif tau_profile_source_inter = tau_profile_source ; ------------------------------------------------------------------------------------ ; Compute age spectra for each individual WAS measurement above 340K. ; ------------------------------------------------------------------------------------ ;bl_g_ind_all = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm,nciv) & bl_g_nh_ind = bl_g_ind_all & bl_g_tr_ind = bl_g_ind_all ;tf_ind_all = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm,nciv) ;for aa = 0, nm-1 do begin ; print,aa ; tii = interpol(indgen(nco2_ext),co2_time_ext,er2_yrfrac[aa]-life_grid_linear/365.) ; print,tii[0],ti[0],er2_yrfrac[aa] ; zii = zi[aa] ; if zii gt nzth-1 then zii = nzth-1 ; if zii lt 6 then zii = 6 ; if zii gt 18 then zii = 18 ; ; for i = 1, nciv-1 do begin ; if er2_was[ci_vars_er2[i],zti[aa]] gt 0 then begin ; for ki = 0, ng-1 do begin ; gr_theta_30yr = g_linear_theta_all[ti_30yr,ki,round(zii)] ; gr_theta_30yr /= total(gr_theta_30yr) ; gr_sub1 = weight_nh * gr_theta_30yr ; gr_sub2 = weight_tr * gr_theta_30yr ; for pp = 0, npeaks_source-1 do begin ; bl_tr_int = interpolate(surf[*,pp,0,0,i],tii) ; for pp2 = 0, npeaks_source-1 do begin ; bl_nh_int = interpolate(surf[*,pp2,0,1,i],tii) ;; bl_nh_int *= scale_nh[aa,i] ; bl_g_nh_ind[pp,pp2,ki,aa,i] = total(gr_sub1 * bl_nh_int) ; bl_g_tr_ind[pp,pp2,ki,aa,i] = total(gr_sub2 * bl_tr_int) ; bl_g_ind_all[pp,pp2,ki,aa,i] = total(gr_sub1 * bl_nh_int) + total(gr_sub2 * bl_tr_int) ; endfor ; endfor ; if aa eq 0 and ki eq 100 then print,i,bl_g_ind_all[0,0,ki,aa,i],bl_g_tr_ind[0,0,ki,aa,i],bl_g_nh_ind[0,0,ki,aa,i],scale_nh[aa,i] ; endfor ; tf_ind_all[*,*,*,aa,i] = er2_was[ci_vars_er2[i],zti[aa]] / bl_g_ind_all[*,*,*,aa,i] ; endif ; endfor ;endfor ;save,bl_g_nh_ind,bl_g_tr_ind,bl_g_ind_all,tf_ind_all,filename=dir+'Aircraft/Missions/seac4rs/bl_tf_ind_'+shapes[shp_i]+'.sav' ;restore,dir+'Aircraft/Missions/seac4rs/bl_tf_ind_'+shapes[shp_i]+'.sav' restore,dir+'Aircraft/Missions/seac4rs/bl_tf_ind.sav' ;co2_g_ind = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm) & co2_g_nh_ind = co2_g_ind & co2_g_ind_diffs = co2_g_ind ;for aa = 0, nm-1 do begin ; if co2_ind[aa] gt 0 then begin ; tii = interpol(indgen(nco2_ext),co2_time_ext,er2_yrfrac[aa]-life_grid_linear/365.) ; zii = zi[aa] ; if zii gt nzth-1 then zii = nzth-1 ; if zii lt 6 then zii = 6 ; if zii gt 18 then zii = 18 ; for ki = 0, ng-1 do begin ; gr_theta_30yr = g_linear_theta_all[ti_30yr,ki,round(zii)] ; gr_theta_30yr /= total(gr_theta_30yr) ; gr_sub1 = weight_nh * gr_theta_30yr ; gr_sub2 = weight_tr * gr_theta_30yr ; for pp = 0, npeaks_source-1 do begin ; co2_tr_int = interpolate(co2_surf[*,pp,0,0],tii) ; for pp2 = 0, npeaks_source-1 do begin ; co2_nh_int = interpolate(co2_surf[*,pp2,0,1],tii) ; co2_g_nh_ind[pp,pp2,ki,aa] = total(gr_sub1 * co2_nh_int) ; co2_g_ind[pp,pp2,ki,aa] = total(gr_sub1 * co2_nh_int) + total(gr_sub2 * co2_tr_int) ; endfor ; endfor ; endfor ; co2_g_ind_diffs[*,*,*,aa] = co2_g_ind[*,*,*,aa] - co2_ind[aa] ; endif ;endfor ;save,co2_g_ind,co2_g_nh_ind,co2_g_ind_diffs,filename=dir+'Aircraft/Missions/seac4rs/co2_diffs_ind_'+shapes[shp_i]+'.sav' ;restore,filename=dir+'Aircraft/Missions/seac4rs/co2_diffs_ind_'+shapes[shp_i]+'.sav' restore,filename=dir+'Aircraft/Missions/seac4rs/co2_diffs_ind.sav' tau_ind = replicate(!values.f_nan,nciv,nm) & diffs_ind = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm,nciv) & ki_ind = replicate(!values.f_nan,nm) & tr_peak_ind = ki_ind diffs_ind_tot = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nm) & sspb_ind = diffs_ind_tot & nh_peak_ind = ki_ind & g_best_ind = fltarr(ngtl,nm) & tfk_best_ind = g_best_ind bl_best_ind = replicate(!values.f_nan,nm,nciv) & bl_tr_best_ind = bl_best_ind & bl_nh_best_ind = bl_best_ind & tf_best_ind = bl_best_ind & co2_best_ind = ki_ind & mean_age_ind = ki_ind mean_age_ind_tr = ki_ind & mean_age_ind_nh = ki_ind & modal_age_ind = ki_ind & modal_age_ind_tr = ki_ind & modal_age_ind_nh = ki_ind & g_best_ind_tr = g_best_ind & g_best_ind_nh = g_best_ind scale_tf = replicate(1.0,nciv,nm) & scale_tf_tau = replicate(1.0,nciv,nm) & tau_ind_init = tau_ind & tr_to_nh_g_ratio = ki_ind & tr_source_frac = ki_ind & tfk_best_ind_init = g_best_ind tf_ind_all_init = tf_ind_all & tf_best_ind_init = bl_best_ind & tf_ind_best_mean_rel = bl_best_ind & tau_ind_best_mean_rel = bl_best_ind & tr_fracs_ind = replicate(!values.f_nan,ng,nm) tf_ind_all_scale = replicate(!values.f_nan,npeaks_source,npeaks_source,ng,nciv) for aa = 0, nm-1 do begin mr_ind = replicate(!values.f_nan,nciv) for i = 1, nciv-1 do mr_ind[i] = er2_was[ci_vars_er2[i],zti[aa]] gd = where(mr_ind[1:-1] gt 0,ngd) if ngd gt 9 then begin ; print,aa zii = zi[aa] if zii gt 18 then zii = 18 if zii lt 6 then zii = 6 ; restore,dir+'Aircraft/Missions/seac4rs/co2_predef_'+shapes[shp_i_opt[zii]]+'.sav' ; restore,dir+'Aircraft/Missions/seac4rs/bl_predef_'+shapes[shp_i_opt[zii]]+'.sav' ; ; restore,dir+'Aircraft/Missions/seac4rs/bl_tf_ind_'+shapes[shp_i_opt[zii]]+'.sav' ; restore,dir+'Aircraft/Missions/seac4rs/co2_diffs_ind_'+shapes[shp_i_opt[zii]]+'.sav' ; for ki = 0, ng-1 do tr_fracs_ind[ki,aa] = total(weight_tr * g_linear_theta_all[ti_30yr,ki,round(zii)]) / total(g_linear_theta_all[ti_30yr,ki,round(zii)]) for i = 1, nciv-1 do begin if mr_ind[i] gt 0 then begin tf_ind_all_scale[*,*,*,i] = tf_ind_all[*,*,*,aa,i] / scale_nh[aa,i] ; temp = tf_ind_all[*,*,*,aa,i] / scale_nh[aa,i] ; chk = where(temp gt 1,nchk) ; if nchk gt 1 then temp[chk] = 1.0 tau_ind[i,aa] = interpolate(tau_profile_source[i,*],zii) endif endfor tau = tau_ind[gd+1,aa] tau_i = interpol(indgen(ngtl),life_grid_linear,tau) for ki = 0, ng-1 do begin ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,round(zii)],tau_i) ; if aa eq 0 and ki eq 64 then print,zii,tfk_i,reform(tf_ind_all[0,0,ki,aa,gd+1]) for pp2 = 0, npeaks_source-1 do for pp = 0, npeaks_source-1 do begin diffs_ind[pp,pp2,ki,aa,gd+1] = (tf_ind_all_scale[pp,pp2,ki,gd+1] - tfk_i)^2 diffs_ind_tot[pp,pp2,ki,aa] = total(diffs_ind[pp,pp2,ki,aa,gd+1]) median_ratio = median(alog(tf_ind_all_scale[pp,pp2,ki,gd+1] / tfk_i)) sspb_ind[pp,pp2,ki,aa] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor if co2_ind[aa] gt 0 then diffs_ind_tot[*,*,ki,aa] = 10. * abs(co2_g_ind_diffs[*,*,ki,aa]) + abs(sspb_ind[*,*,ki,aa]) else diffs_ind_tot[*,*,ki,aa] = abs(sspb_ind[*,*,ki,aa]) endfor ; if aa eq 0 then print,ki_prof[zii[aa]] ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min for ki = 0, ng-1 do begin diff_tot = diffs_ind_tot[*,*,ki,aa] nh_peak_tot = fix(indgen(nsources) / npeaks_source) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[0]] nh_p[ki] = nh_peak_tot[d_tot_si[0]] tr_p[ki] = tr_peak_tot[d_tot_si[0]] diff_min[ki] = d_tot_si[0] ; if aa eq 0 and ki eq 64 then print,diff_tot,d_tot_si[0],diff_val[ki],nh_p[ki],tr_p[ki] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_ind[aa] = ki_min tr_peak_ind[aa] = tr_p[ki_min] nh_peak_ind[aa] = nh_p[ki_min] tfk_best_ind_init[*,aa] = tfk[*,ki_min,round(zii)] tfk_best_ind[*,aa] = tfk_best_ind_init[*,aa] tf_best_ind_init[aa,*] = tf_ind_all_scale[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,*] tfi = interpol(indgen(ngtl),tfk_best_ind[*,aa],tf_best_ind_init[aa,*]) ; tfi = interpol(indgen(ngtl),tfk_best_ind[*,aa],tf_ind_adj[aa,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_ind_init[i,aa] = tau endif endfor tau_ind[*,aa] = tau_ind_init[*,aa] ; tf_best_ind[aa,*] = tf_best_ind_init[aa,*] tf_best_ind[aa,*] = tf_ind_adj[aa,*] ; if aa eq 28 then print,reform(tf_best_ind_init[aa,1:10]),reform(tf_ind_adj[aa,1:10]),tau_ind_init[1:10,aa],co2_ind[aa] ; Check for BL fracs gt 1 and if so redo above diff minimization. chk = where(tf_best_ind[aa,*] ge 1,nchk) if nchk gt 0 then begin ; print,aa,nchk tf = fltarr(nciv) for i = 1, nciv-1 do tf[i] = interpolate(tf_best_source[*,i],zii) scale_tf[chk,aa] = tf[chk]/reform(tf_best_ind[aa,chk]) for ki = 0, ng-1 do begin ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,round(zii)],tau_i) for pp2 = 0, npeaks_source-1 do for pp = 0, npeaks_source-1 do begin tf_ind_all_scale[pp,pp2,ki,*] *= scale_tf[*,aa] diffs_ind[pp,pp2,ki,aa,gd+1] = (tf_ind_all_scale[pp,pp2,ki,gd+1] - tfk_i)^2 diffs_ind_tot[pp,pp2,ki,aa] = total(diffs_ind[pp,pp2,ki,aa,gd+1]) median_ratio = median(alog(tf_ind_all_scale[pp,pp2,ki,gd+1] / tfk_i)) sspb_ind[pp,pp2,ki,aa] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor if co2_ind[aa] gt 0 then diffs_ind_tot[*,*,ki,aa] = 10. * abs(co2_g_ind_diffs[*,*,ki,aa]) + abs(sspb_ind[*,*,ki,aa]) else diffs_ind_tot[*,*,ki,aa] = abs(sspb_ind[*,*,ki,aa]) endfor ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min for ki = 0, ng-1 do begin diff_tot = diffs_ind_tot[*,*,ki,aa] nh_peak_tot = fix(indgen(nsources) / npeaks_source) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[0]] nh_p[ki] = nh_peak_tot[d_tot_si[0]] tr_p[ki] = tr_peak_tot[d_tot_si[0]] diff_min[ki] = d_tot_si[0] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_ind[aa] = ki_min tr_peak_ind[aa] = tr_p[ki_min] nh_peak_ind[aa] = nh_p[ki_min] endif tf_best_ind[aa,*] = tf_ind_all_scale[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,*] ; if aa eq 8 then print,tf_best_ind[aa,3],scale_tf[3,aa] g_best_ind[*,aa] = g_linear_theta_all[*,ki_min,round(zii)] g_best_ind_tr[*,aa] = weight_tr * g_best_ind[*,aa] g_best_ind_nh[*,aa] = weight_nh * g_best_ind[*,aa] ; Adjust the lifetimes to best fit the spectra. bl_best_ind[aa,*] = bl_g_ind_all[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] bl_nh_best_ind[aa,*] = bl_g_nh_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] bl_tr_best_ind[aa,*] = bl_g_tr_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] tfk_best_ind[*,aa] = tfk[*,ki_min,round(zii)] if co2_ind[aa] gt 0 then co2_best_ind[aa] = co2_g_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa] tfi = interpol(indgen(ngtl),tfk_best_ind[*,aa],tf_best_ind[aa,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_ind[i,aa] = tau endif endfor ; tau_ind[*,aa] = tau_ind_init[*,aa] ; if aa eq 28 then print,reform(tf_best_ind[aa,1:10]),tau_ind[1:10,aa] ; Check for anomalous lifetimes and redo the diff minimization if found. for x = 0, 2 do begin tau = fltarr(nciv) for i = 1, nciv-1 do tau[i] = interpolate(tau_profile_source[i,*],zii) ratio = tau_ind[*,aa]/tau chk = where(ratio ge 4. or ratio le 0.25,nchk) if nchk gt 0 then begin tf = fltarr(nciv) for i = 1, nciv-1 do tf[i] = interpolate(tf_best_source[*,i],zii) scale_tf_tau[chk,aa] = tf[chk]/reform(tf_best_ind[aa,chk]) ; low = where(scale_tf_tau[chk,aa] lt 0.9,nlow) ; if nlow gt 0 then scale_tf_tau[chk[low],aa] = 0.9 ; hi = where(scale_tf_tau[chk,aa] gt 1.1,nhi) ; if nhi gt 0 then scale_tf_tau[chk[hi],aa] = 1.1 ; if aa eq 246 then print,scale_tf_tau[chk,aa] for ki = 0, ng-1 do begin ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,round(zii)],tau_i) for pp2 = 0, npeaks_source-1 do for pp = 0, npeaks_source-1 do begin tf_ind_all_scale[pp,pp2,ki,*] *= scale_tf_tau[*,aa] diffs_ind[pp,pp2,ki,aa,gd+1] = (tf_ind_all_scale[pp,pp2,ki,gd+1] - tfk_i)^2 diffs_ind_tot[pp,pp2,ki,aa] = total(diffs_ind[pp,pp2,ki,aa,gd+1]) median_ratio = median(alog(tf_ind_all_scale[pp,pp2,ki,gd+1] / tfk_i)) sspb_ind[pp,pp2,ki,aa] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor if co2_ind[aa] gt 0 then diffs_ind_tot[*,*,ki,aa] = 10. * abs(co2_g_ind_diffs[*,*,ki,aa]) + abs(sspb_ind[*,*,ki,aa]) else diffs_ind_tot[*,*,ki,aa] = abs(sspb_ind[*,*,ki,aa]) endfor ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min for ki = 0, ng-1 do begin diff_tot = diffs_ind_tot[*,*,ki,aa] nh_peak_tot = fix(indgen(nsources) / npeaks_source) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[0]] nh_p[ki] = nh_peak_tot[d_tot_si[0]] tr_p[ki] = tr_peak_tot[d_tot_si[0]] diff_min[ki] = d_tot_si[0] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_ind[aa] = ki_min tr_peak_ind[aa] = tr_p[ki_min] nh_peak_ind[aa] = nh_p[ki_min] tf_best_ind[aa,*] = tf_ind_all_scale[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,*] g_best_ind[*,aa] = g_linear_theta_all[*,ki_min,round(zii)] g_best_ind_tr[*,aa] = weight_tr * g_best_ind[*,aa] g_best_ind_nh[*,aa] = weight_nh * g_best_ind[*,aa] ; Adjust the lifetimes to best fit the spectra. bl_best_ind[aa,*] = bl_g_ind_all[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] bl_nh_best_ind[aa,*] = bl_g_nh_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] bl_tr_best_ind[aa,*] = bl_g_tr_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa,*] tfk_best_ind[*,aa] = tfk[*,ki_min,round(zii)] if co2_ind[aa] gt 0 then co2_best_ind[aa] = co2_g_ind[tr_peak_ind[aa],nh_peak_ind[aa],ki_min,aa] tfi = interpol(indgen(ngtl),tfk_best_ind[*,aa],tf_best_ind[aa,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_ind[i,aa] = tau endif endfor endif endfor ; if aa eq 8 then print,tf_best_ind[aa,3],scale_tf_tau[3,aa] tr_to_nh_g_ratio[aa] = total(g_best_ind_tr[*,aa]) / total(g_best_ind_nh[*,aa]) tr_source_frac[aa] = total(g_best_ind_tr[*,aa]) / total(g_best_ind[*,aa]) mean_age_ind[aa] = total(life_grid_linear*g_best_ind[*,aa]) mean_age_ind_tr[aa] = total(life_grid_linear*g_best_ind_tr[*,aa]) mean_age_ind_nh[aa] = total(life_grid_linear*g_best_ind_nh[*,aa]) mm = max(g_best_ind[*,aa],mi) modal_age_ind[aa] = life_grid_linear[mi] mm = max(g_best_ind_tr[*,aa],mi) modal_age_ind_tr[aa] = life_grid_linear[mi] mm = max(g_best_ind_nh[*,aa],mi) modal_age_ind_nh[aa] = life_grid_linear[mi] endif endfor tf_best_init_ratio = tf_best_ind/tf_ind_init ;bl_nh_scale = 1./tf_best_init_ratio bl_nh_scale = scale_nh for i = 1, nciv-1 do begin ; gd = where(finite(tf_best_init_ratio[*,i]) and scale_nh[*,i] ne 1 and reform(scale_tf_tau[i,*]) ne 1) ; gd = where(scale_nh[*,i] gt 1) ; gd = where(bl_nh_scale[*,i] gt 1) ; bl_nh_scale[gd,i] = 1./(tf_best_init_ratio[gd,i]*(1.-tr_source_frac[gd])) bl_nh_scale[*,i] = scale_nh[*,i] / (reform(scale_tf_tau[i,*]) * reform(scale_tf[i,*])) gd = where(bl_nh_scale[*,i] ne 1) bl_nh_scale[gd,i] /= (1.-tr_source_frac[gd]) gd = where(scale_nh[*,i] lt 1) ; gd = where(bl_nh_scale[*,i] lt 1) bl_nh_scale[gd,i] *= (1.-tr_source_frac[gd]) gd = where(scale_nh[*,i] eq 1) bl_nh_scale[gd,i] = 1. endfor tf_profile_ind = replicate(!values.f_nan,nzth,nciv) & ci_theta_grid_ind_adj = ci_theta_grid_adj for z = 6, 11 do begin zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) for i = 1, nciv-1 do tf_profile_ind[z,i] = mean(tf_best_ind[zi_th,i],/nan) for i = 1, nciv-1 do if finite(tf_profile_ind[z,i]) then ci_theta_grid_ind_adj[z,i] *= (tf_profile_ind[z,i] / tf_best_source[z,i]) endfor ; Smoothing of the average adjusted profiles. ci_theta_grid_ind_adj[8,1] = 10^(0.5 * (alog10(ci_theta_grid_ind_adj[7,1]) + alog10(ci_theta_grid_ind_adj[9,1]))) ci_theta_grid_ind_adj[6,2] *= 1.4 ci_theta_grid_ind_adj[8,2] = 10^(0.75 * alog10(ci_theta_grid_ind_adj[7,2]) + 0.25 * alog10(ci_theta_grid_ind_adj[11,2])) ci_theta_grid_ind_adj[9,2] = 10^(0.5 * alog10(ci_theta_grid_ind_adj[7,2]) + 0.5 * alog10(ci_theta_grid_ind_adj[11,2])) ci_theta_grid_ind_adj[10,2] = 10^(0.25 * alog10(ci_theta_grid_ind_adj[7,2]) + 0.75 * alog10(ci_theta_grid_ind_adj[11,2])) ci_theta_grid_ind_adj[7,3] = 10^(0.5 * (alog10(ci_theta_grid_ind_adj[6,3]) + alog10(ci_theta_grid_ind_adj[8,3]))) ci_theta_grid_ind_adj[13:-1,3] = !values.f_nan ci_theta_grid_ind_adj[16:18,15] *= [0.995,0.98,0.96] ci_theta_grid_ind_adj[16:18,13] *= [0.995,0.98,0.96] ;ci_theta_grid_ind_adj[16:18,12] *= [0.995,0.99,0.975] ;ci_theta_grid_ind_adj[17:18,10] *= [0.95,0.9] for i = 16, 18 do ci_theta_grid_ind_adj[17:18,i] = reform(ci_theta_grid[0,17:18,i]) ci_theta_grid_ind_adj[17:18,14] = reform(ci_theta_grid[0,17:18,14]) for i = 10, 12 do ci_theta_grid_ind_adj[17:18,i] = reform(ci_theta_grid[0,17:18,i]) for i = 10, 20 do ci_theta_grid_ind_adj[6:18,i] = smooth(ci_theta_grid_ind_adj[6:18,i],3);,/edge_truncate) for i = 13, 15, 2 do ci_theta_grid_ind_adj[6:18,i] = smooth(ci_theta_grid_ind_adj[6:18,i],3);,/edge_truncate) ci_theta_grid_ind_adj[6,19] *= 1.00175 ci_theta_grid_ind_adj[7,19] *= 0.9975 ci_theta_grid_ind_adj[8,20] *= 0.999 ; Average lifetime and BL frac profiles based on the individual measurement derived lifetimes and BL fracs. tau_profile_ind = replicate(!values.f_nan,nciv,nzth) & tr_source_profile_ind = tau_profile_ind & nh_source_profile_ind = tau_profile_ind for z = 6, 20 do begin zi_th = where(theta_ind ge theta_grid[z]-dth and theta_ind lt theta_grid[z]+dth,nn) for i = 0, nciv-1 do begin gd = where(finite(tau_ind[i,zi_th]),ngd) if ngd gt 0 then begin tau_profile_ind[i,z] = mean(tau_ind[i,zi_th[gd]]) tf_profile_ind[z,i] = mean(tf_best_ind[zi_th[gd],i]) tf_ind_best_mean_rel[zi_th[gd],i] = tf_best_ind[zi_th[gd],i] / tf_profile_ind[z,i] tau_ind_best_mean_rel[zi_th[gd],i] = tau_ind[zi_th[gd],i] / tau_profile_ind[z,i] endif endfor tr_frac = tr_source_frac[zi_th] gd = where(finite(tr_frac)) tr_frac[gd] /= total(tr_frac[gd]) tr_source_profile_ind[z] = total(tr_frac[gd] * tr_peaks[tr_peak_ind[zi_th[gd]]]) nh_frac = (1. - tr_frac[gd]) nh_frac /= total(nh_frac) nh_source_profile_ind[z] = total(nh_frac * nh_peaks[nh_peak_ind[zi_th[gd]]]) endfor tau_profile_ind[1,11] = !values.f_nan tf_profile_ind[11,1] = !values.f_nan tau_profile_ind[3,13:14] *= [0.8,0.9] tf_profile_ind[12:20,2] = tf_profile_ind[11,2] tf_profile_ind[14:20,3] = (0.99-findgen(7)*0.01) * tf_profile_ind[13,3] for i = 2, 3 do begin chk = where(~finite(tau_profile_ind[i,6:20])) tau_profile_ind[i,chk+6] = 0.5*(tau_profile_ind[i,chk+5]+tau_profile_ind[i,chk+7]) tf_profile_ind[chk+6,i] = 0.5*(tf_profile_ind[chk+5,i]+tf_profile_ind[chk+7,i]) endfor for i = 2, nciv-1 do begin tau_profile_ind[i,6:20] = smooth(tau_profile_ind[i,6:20],3,/edge_truncate) tf_profile_ind[6:20,i] = smooth(tf_profile_ind[6:20,i],3,/edge_truncate) endfor ;print,reform(tf_ind_all[0,0,100,246,1:-1]/tf_ind_all_init[0,0,100,246,1:-1]) ;save,ci_theta_grid_ind_adj,filename=dir+'Aircraft/Missions/seac4rs/theta_avg_profiles.sav' restore,filename=dir+'Aircraft/Missions/seac4rs/theta_avg_profiles.sav' ; Redo with adjusted BL fracs ; -------------------------------------- Age spectra for theta grid source regions. --------------------------------------------------------------------------------------------------- doit = 0 ; *********************************** if doit then begin for z = 6, 18 do begin gd = where(finite(reform(ci_theta_grid_ind_adj[z,1:-1])),ngd) if z ge 7 then begin tau_profile_source[*,z] = replicate(!values.f_nan,nciv) tau_profile_source[gd+1,z] = tau_profile_source[gd+1,z-1] endif tau = tau_profile_source[gd+1,z] tau_i = interpol(indgen(ngtl),life_grid_linear,tau) for ki = 0, ng-1 do begin for i = 1, nciv-1 do begin tf_surf[*,*,*,*,ki,z,i] = ci_theta_grid_ind_adj[z,i]/bl_g_all[*,*,*,*,ki,z,i] temp = tf_surf[*,*,*,*,ki,z,i] chk = where(temp gt 1,nchk) if nchk gt 1 then temp[chk] = 1.0 tf_surf[*,*,*,*,ki,z,i] = temp endfor ; Find diffs from mu-tau relationships. tfk_i = interpolate(tfk[*,ki,z],tau_i) for wi2 = 0, nwidths_source-1 do for pp2 = 0, npeaks_source-1 do for wi = 0, nwidths_source-1 do for pp = 0, npeaks_source-1 do begin diffs_source[pp,wi,pp2,wi2,ki,z,gd+1] = (tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] - tfk_i)^2 diffs_source_tot[pp,wi,pp2,wi2,ki,z] = total(diffs_source[pp,wi,pp2,wi2,ki,z,gd+1]) median_ratio = median(alog(tf_surf[pp,wi,pp2,wi2,ki,z,gd+1] / tfk_i)) sspb_source[pp,wi,pp2,wi2,ki,z] = 1e2 * signum(median_ratio) * (exp(abs(median_ratio)) - 1.0) endfor diff_tot_source[*,*,*,*,ki,z] = 10. * abs(co2_g_all_diffs[*,*,*,*,ki,z]) + abs(sspb_source[*,*,*,*,ki,z]) si_co2_ki[*,ki,z] = sort(abs(co2_g_all_diffs[*,*,*,*,ki,z])) si_diffs_ki[*,ki,z] = sort(diffs_source_tot[*,*,*,*,ki,z]) si_sspb_ki[*,ki,z] = sort(sspb_source[*,*,*,*,ki,z]) sspb_ki_med[ki,z] = median(sspb_source[*,*,*,*,ki,z]) co2_temp = co2_g_all_diffs[*,*,*,*,ki,z] diffs_temp = diffs_source_tot[*,*,*,*,ki,z] sspb_temp = sspb_source[*,*,*,*,ki,z] ; if z eq 6 then print,co2_temp[si_co2_ki[0:5,ki,z]],sspb_temp[si_sspb_ki[0:5,ki,z]],sspb_ki_med[ki,z] endfor ; Check for best fits to tracers and CO2 for ki values. nsources = npeaks_source^2*nwidths_source^2 diff_min = replicate(!values.f_nan,ng) & diff_val = diff_min & nh_w = diff_min & nh_p = diff_min & tr_w = diff_min & tr_p = diff_min if z le 10 then kd = 3 else kd = 2 if z eq 6 then ki_max = ng-1 else ki_max = ki_prof[z-1] - kd for ki = 0, ki_max do begin diff_tot = diff_tot_source[*,*,*,*,ki,z] nh_width_tot = fix(indgen(nsources)/(npeaks_source^2*nwidths_source)) nh_peak_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source) / (npeaks_source*nwidths_source)) tr_width_tot = fix((indgen(nsources) - nh_width_tot*npeaks_source^2*nwidths_source - nh_peak_tot*npeaks_source*nwidths_source) / (npeaks_source)) tr_peak_tot = fix(indgen(nsources) mod npeaks_source) si0 = 0 if z ge 7 then begin si = sort(diff_tot) trp = fix(tr_peak[z-1]) si1 = si[0] chk = where(tr_peak_tot ge trp-1 and tr_peak_tot le trp) d_tot_si = sort(diff_tot[chk]) si0 = d_tot_si[0] diff_val[ki] = diff_tot[chk[si0]] nh_w[ki] = nh_width_tot[chk[si0]] nh_p[ki] = nh_peak_tot[chk[si0]] tr_w[ki] = tr_width_tot[chk[si0]] tr_p[ki] = tr_peak_tot[chk[si0]] endif else begin d_tot_si = sort(diff_tot) diff_val[ki] = diff_tot[d_tot_si[si0]] nh_w[ki] = nh_width_tot[d_tot_si[si0]] nh_p[ki] = nh_peak_tot[d_tot_si[si0]] tr_w[ki] = tr_width_tot[d_tot_si[si0]] tr_p[ki] = tr_peak_tot[d_tot_si[si0]] endelse diff_min[ki] = d_tot_si[si0] endfor dsort = sort(diff_val) dmin = min(diff_val,ki_min) ki_prof[z] = ki_min tr_peak[z] = tr_p[ki_min] tr_width[z] = tr_w[ki_min] nh_peak[z] = nh_p[ki_min] nh_width[z] = nh_w[ki_min] g_best_source[*,z] = g_linear_theta_all[*,ki_min,z] g_best_source_tr[*,z] = weight_tr * g_best_source[*,z] g_best_source_nh[*,z] = weight_nh * g_best_source[*,z] tr_best_source_frac[z] = total(g_best_source_tr[*,z]) / total(g_best_source[*,z]) ; Adjust the lifetimes to best fit the spectra. bl_best_source[z,*] = bl_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_nh_best_source[z,*] = bl_g_nh[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] bl_tr_best_source[z,*] = bl_g_tr[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] tfk_best_source[*,z] = tfk[*,ki_min,z] tf_best_source[z,*] = tf_surf[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z,*] co2_best_source[z] = co2_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_min,z] tfi = interpol(indgen(ngtl),tfk_best_source[*,z],tf_best_source[z,*]) for i = 1, nciv-1 do begin if finite(tfi[i]) then begin tau = interpolate(life_grid_linear,tfi[i]) tau_profile_source[i,z] = tau endif endfor mean_age_source[z] = total(life_grid_linear*g_best_source[*,z]) mean_age_source_tr[z] = total(life_grid_linear*g_best_source_tr[*,z]) mean_age_source_nh[z] = total(life_grid_linear*g_best_source_nh[*,z]) mm = max(g_best_source[*,z],mi) modal_age_source[z] = life_grid_linear[mi] mm = max(g_best_source_tr[*,z],mi) modal_age_source_tr[z] = life_grid_linear[mi] mm = max(g_best_source_nh[*,z],mi) modal_age_source_nh[z] = life_grid_linear[mi] endfor ; print,ki_prof[5:18],tr_peak[5:18],nh_peak[5:18] endif zii = zi zii[where(zii gt 18)] = 18 zii[where(zii lt 6)] = 6 bl_nh_ind_scaled = bl_nh_scale for i = 1, nciv-1 do begin gd = where(bl_nh_scale[*,i] ne 1) bl_nh_ind_scaled[gd,i] = bl_nh_scale[gd,i] * bl_nh_best_source[zii[gd],i] / (1.0-tr_best_source_frac[zii[gd]]) gd = where(bl_nh_scale[*,i] eq 1) bl_nh_ind_scaled[gd,i] = !values.f_nan endfor ; Read in age spectra from modeling studies. data = read_ascii(dir+'Models/CLaMS/Hauck_spectrum_JJA_370K_56N.txt',data_start=1) age_1 = reform(data.field1[0,*]) freq_1 = reform(data.field1[1,*]) freq_1 /= total(freq_1) data = read_ascii(dir+'Models/CLaMS/Hauck_spectrum_SON_370K_56N.txt',data_start=1) age_2 = reform(data.field1[0,*]) freq_2 = reform(data.field1[1,*]) freq_2 /= total(freq_2) data = read_ascii(dir+'Models/CLaMS/Hauck_age_spectrum_JJA_70hPa_55N.txt',data_start=1) age_3 = reform(data.field1[0,*]) freq_3 = reform(data.field1[1,*]) freq_3 /= total(freq_3) data = read_ascii(dir+'Models/CLaMS/Hauck_age_spectrum_SON_70hPa_55N.txt',data_start=1) age_4 = reform(data.field1[0,*]) freq_4 = reform(data.field1[1,*]) freq_4 /= total(freq_4) data = read_ascii(dir+'Models/CLaMS/Ploeger_age_spectrum_JJA_400K_60N.txt',data_start=1) age_5 = reform(data.field1[0,*]) freq_5 = reform(data.field1[1,*]) freq_5 /= total(freq_5) data = read_ascii(dir+'Models/CLaMS/Ploeger_age_spectrum_SON_400K_60N.txt',data_start=1) age_6 = reform(data.field1[0,*]) freq_6 = reform(data.field1[1,*]) freq_6 /= total(freq_6) ; Save the data for publication. ;gatt = hash('Title','Data used in Ray et al., Age spectra and other transport diagnostics in the North American monsoon upper troposphere and lower stratosphere from SEAC4RS in situ measurements') ;dim_names = ['ny_surface','nt_surface','ntracers','nz_theta_grid'] ;atts = ['Surface latitude grid','Surface times','Theta grid','Surface data interpolated in latitude','CO2 Surface data from CarbonTracker interpolated in latitude', $ ; 'Theta average trace gas profiles','Theta average CO2 profiles'] ;vars = hash('Lat_Surface',hash('value',lat_surf_interp,'dim_names',dim_names[0],'dim_sizes',ny_surf_interp,'attributes',hash('description',atts[0],'units','Degrees')), $ ; 'Time_Surface',hash('value',co2_time_ext,'dim_names',dim_names[1],'dim_sizes',nco2_ext,'attributes',hash('description',atts[1],'units','Years')), $ ; 'Theta_Grid',hash('value',theta_grid,'dim_names',dim_names[3],'dim_sizes',nzth,'attributes',hash('description',atts[2],'units','K')), $ ; 'Surface_Data',hash('value',surface_data_interp[*,*,1:-1],'dim_names',[dim_names[1],dim_names[0],dim_names[2]],'dim_sizes',[nco2_ext,ny_surf_interp,nciv-1],'attributes', $ ; hash('description',atts[3],'units','mixing ratio')), $ ; 'Surface_Data_CO2',hash('value',co2_surface_interp,'dim_names',[dim_names[1],dim_names[0]],'dim_sizes',[nco2_ext,ny_surf_interp],'attributes',hash('description',atts[4],'units','ppm')), $ ; 'Tracer_Theta_Avgs',hash('value',reform(ci_theta_grid[0,*,1:-1]),'dim_names',[dim_names[3],dim_names[2]],'dim_sizes',[nzth,nciv-1],'attributes',hash('description',atts[5],'units','mixing ratio')), $ ; 'CO2_Theta_Avgs',hash('value',reform(er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')]),'dim_names',[dim_names[3]],'dim_sizes',[nzth],'attributes',hash('description',atts[6], $ ; 'units','ppm'))) ;ncdf_put,'/Users/eray/Work/Papers/SEAC4RS Ozone/SEAC4RS_UTLS_transport_pub_data.nc',variables=vars,/new,gatt=gatt ; ;openw,1,'/Users/eray/Work/Papers/SEAC4RS Ozone/SEAC4RS_UTLS_transport_pub_tracers.txt' ;printf,1,'Trace gases used in Ray et al., Age spectra and other transport diagnostics in the North American monsoon upper troposphere and lower stratosphere from SEAC4RS in situ measurements' ;for i = 1, nciv-1 do printf,1,ci_vars_dc8[i] ;close,1 ;id = ncdf_create('/Users/eray/Work/Papers/SEAC4RS Ozone/SEAC4RS_UTLS_transport_pub_data2.nc',/clobber,/netcdf4_format) ;dimid = intarr(4) ;dimid[0] = ncdf_dimdef(id,'ny_surface',ny_surf_interp) ;dimid[1] = ncdf_dimdef(id,'nt_surface',nco2_ext) ;dimid[2] = ncdf_dimdef(id,'nz_theta_grid',nzth) ;dimid[3] = ncdf_dimdef(id,'ntracers',nciv-1) ;id0 = ncdf_vardef(id,'Lat_Surface',[dimid[0]],/float) ;id1 = ncdf_vardef(id,'Time_Surface',[dimid[1]],/float) ;id2 = ncdf_vardef(id,'Theta_Grid',[dimid[2]],/float) ;id4 = ncdf_vardef(id,'Surface_Data',[dimid[1],dimid[0],dimid[3]],/float) ;id5 = ncdf_vardef(id,'Surface_Data_CO2',[dimid[1],dimid[0]],/float) ;id6 = ncdf_vardef(id,'Tracer_Theta_Avgs',[dimid[2],dimid[3]],/float) ;ncdf_attput,id,/global,'Title', $ ; 'Data used in Ray et al., Age spectra and other transport diagnostics in the North American monsoon upper troposphere and lower stratosphere from SEAC4RS in situ measurements' ;ncdf_control,id,/endef ;ncdf_varput,id,id0,lat_surf_interp ;ncdf_varput,id,id1,co2_time_ext ;ncdf_varput,id,id2,theta_grid ;ncdf_varput,id,id4,reform(surface_data_interp[*,*,1:-1]) ;ncdf_varput,id,id5,co2_surface_interp ;ncdf_varput,id,id6,reform(ci_theta_grid[0,*,1:-1]) ;ncdf_close,id colors = ['purple','blue','sky blue','lime green','green','orange','orange red','black','gray','sienna', $ 'purple','blue','sky blue','lime green','green','orange','orange red','black','gray','sienna'] cgloadct,3,rgb_table=rgb_3 cgloadct,30,rgb_table=rgb_30 cgloadct,39,rgb_table=rgb_39 ;ord_min = 1000 & ord_max = 3 ;logaxis_ticks,ord_min,ord_max,tickvals ;yticks = n_elements(tickvals) ;yticknames = strcompress(string(tickvals)) ;ytickv = tickvals ;yticknames[1:6] = replicate('',6) & yticknames[8:9] = replicate('',2) & yticknames[11:16] = replicate('',6) & yticknames[18] = '' & yticknames[20:25] = replicate('',6) ;yti = [0,7,10,17,19,26] ; ;p = plot(indgen(2),/nodata,xrange=[-90,90],xminor=0,xtitle='Latitude',ytitle='Pressure (hPa)',yrange=[ord_min+20,ord_max],/ystyle,/ylog,ymajor=6,yminor=0,ytickname=reverse(yticknames[yti]), $ ; ytickv=ytickv[yti],xticklen=0.02,yticklen=0.02,xtickinterval=30,axis_style=1,font_size=11,dimensions=[800,500]) ;p = plot(tpzm['lat','value'],tpause,thick=2,color='orange',/overplot) ;;c = contour(theta_mean[*,19:-1],tpzm['lat','value'],pssr[19:-1],c_thick=2,c_value=[380,380],c_color=['orange','orange'],c_linestyle=1,c_label_show=0,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_zm.png' ord_min = 1000 & ord_max = 30 logaxis_ticks,ord_min,ord_max,tickvals yticks = n_elements(tickvals) yticknames = strcompress(string(tickvals)) ytickv = tickvals yticknames[1:6] = replicate('',6) & yticknames[8:9] = replicate('',2) & yticknames[11:16] = replicate('',6) yti = [0,7,10,17] ;p = plot(indgen(2),/nodata,xrange=[-90,90],xminor=0,xtitle='Latitude',ytitle='Pressure (hPa)',yrange=[ord_min+20,ord_max],/ystyle,/ylog,ymajor=6,yminor=0,ytickname=reverse(yticknames[yti]), $ ; ytickv=ytickv[yti],xticklen=0.02,yticklen=0.02,xtickinterval=30,axis_style=1,font_size=11,dimensions=[800,500]) ;;p = plot(tpzm['lat','value'],tpause,thick=2,color='orange',/overplot) ;p = plot(tpzm['lat','value'],tp_m,color='orange',thick=2,/overplot) ;;poly = polygon([tpzm['lat','value',220:-1],tpzm['lat','value',-1],tpzm['lat','value',220]],[tpause[220:-1],100,100],/data,/fill_background,fill_color='red',fill_transparency=80,color='white') ;;poly = polygon([tpzm['lat','value',0:140],tpzm['lat','value',140],tpzm['lat','value',0]],[tpause[0:140],100,100],/data,/fill_background,fill_color='red',fill_transparency=80,color='white') ;c = contour(theta_mean[*,19:-1],tpzm['lat','value'],pssr[19:-1],c_thick=2,c_value=[350,360,380,400,450,500],c_color=replicate('orange',6),c_linestyle=1,c_label_show=0,/overplot) ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_zm2.png' ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_zm3.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_zm4.png' ;tii = indgen(10)+500 ;map = MAP('Equirectangular',linestyle='dotted',title=title,limit=[10,230,53,290],label_position=0,/box_axes,margin=0.11,dimensions=[800,500]) ;m1 = MAPCONTINENTS(/countries,fill_color='beige') ;m1 = MAPCONTINENTS(/usa) ;zi = where(altp_er2 ge 12) ;p = plot(lon_er2[zi],lat_er2[zi],color='blue',symbol='o',sym_size=0.25,linestyle=6,/overplot) ;zi = where(altp_er2[gd3] ge 12) ;p = plot(lon_er2[gd3[zi[tii]]],lat_er2[gd3[zi[tii]]],color='orange',symbol='o',sym_size=0.25,linestyle=6,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_map.png' ;p = plot(indgen(2),/nodata,xrange=[-10,65],yrange=[0,0.1],xtitle='Latitude',ytitle='Scaled Frequency',title='Source Latitude Distributions',dimensions=[600,400],margin=[0.11,0.1,0.03,0.08], $ ; font_size=11) ;offset = 0 ;for z = 6, 18 do begin ; tr_tot = total(g_best_source_tr[*,z]) ; if z ge 7 then begin ; if tr_peak[z] eq tr_peak[z-1] then offset -= 0.4 else offset = 0 ; endif ;; p = plot(tr_source_lats+offset,tr_tot*tr_gauss[*,tr_peak[z],tr_width[z]],color='lime green',thick=4,linestyle=1,/overplot) ; p = plot(tr_source_lats+offset,tr_tot*tr_gauss[*,tr_peak[z],tr_width[z]],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;endfor ;tr_tot = total(g_best_source_tr[*,10]) ;p = plot(tr_source_lats,tr_tot*tr_gauss[*,tr_peak[10],tr_width[10]],color='deep sky blue',thick=2,/overplot) ;offset = 0 ;for z = 6, 14 do begin ; nh_tot = total(g_best_source_nh[*,z]) ; if z ge 7 then begin ; if tr_peak[z] eq tr_peak[z-1] then offset -= 0.4 else offset = 0 ; endif ;; p = plot(nh_source_lats+offset,nh_tot*nh_gauss[*,nh_peak[z],nh_width[z]],color='medium orchid',thick=4,linestyle=1,/overplot) ; p = plot(nh_source_lats+offset,nh_tot*nh_gauss[*,nh_peak[z],nh_width[z]],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ; if z eq 10 then p = plot(nh_source_lats+offset,nh_tot*nh_gauss[*,nh_peak[z],nh_width[z]],color='deep sky blue',thick=2,/overplot) ;endfor ;tickvals = findgen(7)*20+theta_grid[6]-5 ;tickname = strmid(strcompress(string(tickvals),/r),0,3) ;cb = COLORBAR(POSITION=[0.9,0.3,0.92,0.75],tickname=tickname,tickvalues=1.5*(tickvals-theta_grid[6]+5),RGB_TABLE=rgb_3[70:-1,*],BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir, $ ; font_size=10,title='Theta (K)') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_source_lat_dists.png' z = 14 ;p = plot(indgen(2),/nodata,xrange=[-20,60],yrange=[0,0.1],xtitle='Latitude',ytitle='Scaled Frequency',title='Source Latitude Distributions',dimensions=[650,400],margin=[0.11,0.1,0.12,0.08], $ ; font_size=11) ;nh_tot = total(g_best_source_nh[*,z]) ;p = plot(nh_source_lats,nh_tot*nh_gauss[*,nh_peak[z],nh_width[z]],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;tr_tot = total(g_best_source_tr[*,z]) ;p = plot(tr_source_lats,tr_tot*tr_gauss[*,tr_peak[z],tr_width[z]],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;for y = 0, nn-1 do begin ; nh_tot = total(g_best_ind_nh[*,zi_th[y]]) ; if nh_tot gt 0 and co2_ind[zi_th[y]] gt 0 then p = plot(nh_source_lats,nh_tot*nh_gauss[*,nh_peak_ind[zi_th[y]],0],color='sky blue',thick=1,/overplot) ; tr_tot = total(g_best_ind_tr[*,zi_th[y]]) ; if tr_tot gt 0 and co2_ind[zi_th[y]] gt 0 then p = plot(tr_source_lats,tr_tot*tr_gauss[*,tr_peak_ind[zi_th[y]],0],color='orange',thick=1,/overplot) ;endfor ;p = plot(indgen(2),/nodata,xrange=[387,400],yrange=[290,500],xtitle='ppmv',ytitle='Theta (K)',title='SEAC$^4$RS CO$_2$',dimensions=[700,500],margin=0.09,font_size=11) ;p = plot(dc8_was['CO2'],dc8_was['THETA'],color='sky blue',symbol='o',sym_size=0.5,/sym_filled,linestyle=6,/overplot) ;p = plot(er2_was['CO2_HUPCRS'],er2_was['THETA'],color='blue',symbol='o',sym_size=0.5,/sym_filled,linestyle=6,/overplot) ;p = errorplot(er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')],theta_grid,er2_was_theta_grid[2,*,where(varnames_er2_was eq 'CO2_HUPCRS')],replicate(0,nzth),color='purple', $ ; thick=3,symbol='o',/sym_filled,errorbar_capsize=0,/overplot) ;z = 6 ;zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;gd = where(finite(tf_ind_init[zi_th,3])) ;si = sort(tf_ind_init[zi_th[gd],3]) ;s = symbol(co2_ind[zi_th[gd[si[-6]]]],theta_ind[zi_th[gd[si[-6]]]],sym_color='medium orchid','o',/sym_filled,/data) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_theta_profile.png' colors3 = ['plum','blue','sky blue','green','orange','red','brown','tan'] ;p = plot(indgen(2),/nodata,xrange=[389.5,397],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='SEAC$^4$RS CO$_2$',dimensions=[500,500],font_size=11,margin=[0.13,0.09,0.05,0.08]) ;;p = errorplot(er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')],theta_grid,er2_was_theta_grid[2,*,where(varnames_er2_was eq 'CO2_HUPCRS')],replicate(0,nzth),color='purple', $ ;; thick=3,symbol='o',/sym_filled,errorbar_capsize=0,/overplot) ;;for b = 7, 14 do begin ;; prof = replicate(!values.f_nan,nzth) ;; for z = 6, 18 do prof[z] = co2_g[b,diff_min_i[z],z] ;; p = plot(prof,theta_grid,color=colors[b-6],symbol='o',/sym_filled,thick=2,/overplot) ;; t = text(0.8,0.81-(b-7.)*0.04,strcompress(string(fix(lat_ct[b]-5)),/r)+'-'+strcompress(string(fix(lat_ct[b]+5)),/r),color=colors[b-6],font_size=11,/norm) ;;endfor ;;prof = replicate(!values.f_nan,nzth) ;for pp = 1, npeaks_source-2 do for pp2 = 0, npeaks_source-3 do begin ;; for z = 6, 18 do prof[z] = co2_g_all[pp,0,pp2,0,ki_prof[z],z] ;; p = plot(prof,theta_grid,color=colors3[pp-1],thick=2,transparency=pp2*10,/overplot) ;;; for z = 6, 18 do prof[z] = co2_g_all[pp,0,5,0,ki_prof[z],z] ;;; p = plot(prof,theta_grid,color=colors[pp+1],thick=2,linestyle=2,/overplot) ;;; for z = 6, 18 do prof[z] = co2_g_all[pp,0,3,0,ki_prof[z],z] ;;; p = plot(prof,theta_grid,color=colors[pp+1],thick=2,linestyle=1,/overplot) ;;; for z = 6, 18 do prof[z] = co2_g_all[pp,0,2,0,ki_prof[z],z] ;;; p = plot(prof,theta_grid,color=colors[pp+1],thick=2,linestyle=3,/overplot) ; if pp2 eq 1 then t = text(0.27,0.41-(pp)*0.04,strcompress(string(fix(tr_peaks[pp])),/r)+'$\deg$N',color=colors3[pp-1],font_size=11,alignment=0.5,/norm) ; if pp eq 1 then t = text(0.43,0.41-(pp2+1)*0.04,strcompress(string(fix(nh_peaks[pp2])),/r)+'$\deg$N',color=colors3[pp2],font_size=11,alignment=0.5,/norm) ;endfor ;p = plot(er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')],theta_grid,color='purple',thick=3,symbol='o',/sym_filled,/overplot) ;;p = plot(co2_g_best,theta_grid,color='medium_orchid',symbol='o',/sym_filled,thick=2,/overplot) ;;for z = 6, 18 do prof[z] = co2_g_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_prof[z],z] ;;p = plot(prof,theta_grid,color='medium_orchid',symbol='o',/sym_filled,thick=2,/overplot) ;for z = 6, 18 do begin ; tr_tot = total(g_best_source_tr[*,z]) ; s = symbol(co2_best_source[z],theta_grid[z],'s',sym_color=colors3[tr_peak[z]-1],/sym_filled,sym_size=tr_tot*2.5,/data) ;endfor ;for z = 6, 14 do begin ; nh_tot = total(g_best_source_nh[*,z]) ; s = symbol(co2_best_source[z],theta_grid[z],'o',sym_color=colors3[nh_peak[z]],sym_size=nh_tot*5,sym_thick=2,/data) ;endfor ;;;;for s = 0, 9, 3 do p = plot(co2_g_best_seas[*,s],theta_grid,color='pink',thick=2,/overplot) ;t = text(0.35,0.54,'CarbonTracker',font_size=11,alignment=0.5,/norm) ;t = text(0.27,0.5,'Tropical',font_size=11,alignment=0.5,/norm) ;t = text(0.27,0.46,'Lat Peak',font_size=11,alignment=0.5,/norm) ;t = text(0.27,0.42,'(squares)',font_size=11,alignment=0.5,/norm) ;t = text(0.43,0.5,'NH',font_size=11,alignment=0.5,/norm) ;t = text(0.43,0.46,'Lat Peak',font_size=11,alignment=0.5,/norm) ;t = text(0.43,0.42,'(circles)',font_size=11,alignment=0.5,/norm) ;;t = text(0.15,0.2,'SEAC$^4$RS Measured',color='purple',font_size=11,/norm) ;;t = text(0.17,0.35,'LEF',color=colors[5],font_size=11,/norm) ;;t = text(0.17,0.32,'SGP',color=colors[4],font_size=11,/norm) ;;t = text(0.17,0.29,'KEY',color=colors[3],font_size=11,/norm) ;;t = text(0.17,0.26,'MLO',color=colors[2],font_size=11,/norm) ;;t = text(0.17,0.23,'SMO',color=colors[1],font_size=11,/norm) ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_gr_predef_site_profiles.png' ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_gr_source_profiles.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_gr_source_profiles1.png' ;p = plot(indgen(2),/nodata,xrange=[389.5,397],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='CO$_2$ colored by peak tropical latitudes',dimensions=[500,500],font_size=11, $ ; margin=[0.13,0.09,0.05,0.08]) ;for pp = 1, npeaks_source-2 do begin ; gd = where(tr_peak_ind eq pp,ngd) ; if ngd gt 0 then begin ; p = plot(co2_ind[gd],theta_ind[gd],symbol='o',color=colors3[pp-1],linestyle=6,/sym_filled,/overplot) ; endif ;endfor ;p = plot(indgen(2),/nodata,xrange=[389.5,397],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='CO$_2$ colored by peak NH latitudes',dimensions=[500,500],font_size=11, $ ; margin=[0.13,0.09,0.05,0.08]) ;for pp = 0, npeaks_source-1 do begin ; gd = where(nh_peak_ind eq pp,ngd) ; if ngd gt 0 then begin ; p = plot(co2_ind[gd],theta_ind[gd],symbol='o',color=colors3[pp],linestyle=6,/sym_filled,/overplot) ; endif ;endfor ;p = plot(indgen(2),/nodata,xrange=[389.5,397],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='CO$_2$ Colored by BL Source Peak Latitudes',dimensions=[500,500],font_size=11, $ ; margin=[0.13,0.09,0.03,0.08]) ;for pp = 0, npeaks_source-1 do begin ; gd = where(tr_peak_ind eq pp,ngd) ; if ngd gt 0 then for aa = 0, ngd-1 do s = symbol(co2_ind[gd[aa]],theta_ind[gd[aa]],'s',sym_color=colors3[pp],/sym_filled,sym_size=tr_source_frac[gd[aa]]*1.5,sym_transparency=50,/data) ; if fix(tr_peaks[pp]) lt 0 then hem = 'S' ; if fix(tr_peaks[pp]) gt 0 then hem = 'N' ; if fix(tr_peaks[pp]) eq 0 then hem = '' ; t = text(0.24,0.45-(pp+1)*0.04,strcompress(string(fix(abs(tr_peaks[pp]))),/r)+'$\deg$'+hem,color=colors3[pp],font_size=11,alignment=0.5,/norm) ;endfor ;for pp = 0, npeaks_source-1 do begin ; gd = where(nh_peak_ind eq pp,ngd) ; if ngd gt 0 then for aa = 0, ngd-1 do if tr_source_frac[gd[aa]] lt 0.8 then s = symbol(co2_ind[gd[aa]],theta_ind[gd[aa]],'o',sym_color=colors3[pp],sym_size=(1.-tr_source_frac[gd[aa]])*3, $ ; sym_thick=2,/data) ; t = text(0.4,0.45-(pp+1)*0.04,strcompress(string(fix(nh_peaks[pp])),/r)+'$\deg$N',color=colors3[pp],font_size=11,alignment=0.5,/norm) ;endfor ;t = text(0.24,0.54,'Tropical',font_size=11,alignment=0.5,/norm) ;t = text(0.24,0.5,'Lat Peak',font_size=11,alignment=0.5,/norm) ;t = text(0.24,0.46,'(squares)',font_size=11,alignment=0.5,/norm) ;t = text(0.4,0.54,'NH',font_size=11,alignment=0.5,/norm) ;t = text(0.4,0.5,'Lat Peak',font_size=11,alignment=0.5,/norm) ;t = text(0.4,0.46,'(circles)',font_size=11,alignment=0.5,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_gr_source_ind.png' ;p = plot(indgen(2),/nodata,xrange=[389.5,397.5],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='SEAC$^4$RS CO$_2$',dimensions=[700,500],margin=0.09,font_size=11) ;prof = replicate(!values.f_nan,nzth) ;for pp = 0, npeaks_source-1 do for pp2 = 0, npeaks_source-1 do for wi2 = 0, nwidths_source-1 do for wi = 0, nwidths_source-1 do begin ; for z = 6, 18 do prof[z] = co2_clams_all[pp,wi,pp2,wi2,z] ; p = plot(prof,theta_grid,color=colors[pp+1],thick=2,transparency=(5-pp2)*20,/overplot) ;endfor ;p = plot(er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')],theta_grid,color='purple',thick=3,symbol='o',/sym_filled,errorbar_capsize=0,/overplot) ;p = plot(indgen(2),/nodata,xrange=[-3,3],yrange=[345,485],xtitle='ppmv',ytitle='Theta (K)',title='SEAC$^4$RS CO$_2$ Differences',dimensions=[550,550],margin=[0.11,0.07,0.02,0.07],font_size=11) ;p = plot([0,0],[340,490],linestyle=1,/overplot) ;for pp = tr_peak[6]-1, tr_peak[6] do p = plot(co2_g_all_diffs[pp,0,*,0,ki_prof[6],6],replicate(theta_grid[6],npeaks_source),symbol='o',/sym_filled,color=colors3[pp-1], $ ; linestyle=6,/overplot) ;for z = 7, 18 do for pp = tr_peak[z-1]-1, tr_peak[z-1] do p = plot(co2_g_all_diffs[pp,0,*,0,ki_prof[z],z],replicate(theta_grid[z],npeaks_source),symbol='o',/sym_filled,color=colors3[pp-1], $ ; linestyle=6,/overplot) ;for z = 6, 18 do s = symbol(co2_g_all_diffs[tr_peak[z],0,nh_peak[z],0,ki_prof[z],z],theta_grid[z],'s',/sym_filled,sym_color=colors3[tr_peak[z]-1],sym_size=1.5,/data) z = 9 th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;p = plot(indgen(2),/nodata,yrange=[0,0.4],dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([ki_prof[z]-10,ki_prof[z]+10],[0,0],linestyle=1,/overplot) ;for pp = tr_peak[z-1]-1, tr_peak[z-1] do for ki = ki_prof[z]-10, ki_prof[z]+10 do p = plot(replicate(ki,npeaks_source),reform(abs(co2_g_all_diffs[pp,0,*,0,ki,z])),symbol='o',/sym_filled, $ ; sym_size=1,color=colors3[pp-1],linestyle=6,/overplot) ;s = symbol(ki_prof[z],abs(co2_g_all_diffs[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_prof[z],z]),symbol='s',/sym_filled,sym_size=1.5,sym_color='blue',/data) ;p = plot(indgen(2),/nodata,xrange=[Kdiff[ki_prof[z]-30],1e3],/xlog,yrange=[0,32],xtitle='Diffusion coefficient (K,m$^2$/s)',ytitle='Differences', $ ; title='Measured-Idealized BL Frac and CO$_2$ Differences vs. Diffusion Coefficient '+th1+'K',dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([ki_prof[z]-30,ng-1],[0,0],linestyle=1,/overplot) ;min_sspb = fltarr(ng) & max_sspb = min_sspb & min_co2 = min_sspb & max_co2 = min_sspb & min_tot = min_sspb & max_tot = min_sspb ;if z ge 7 then ppi = indgen(2)+tr_peak[z-1]-1 else ppi = indgen(npeaks_source) ;for ki = 0, ng-1 do begin ; min_sspb[ki] = min(abs(sspb_source[ppi,0,*,0,ki,z])) ; max_sspb[ki] = max(abs(sspb_source[ppi,0,*,0,ki,z])) ; min_co2[ki] = min(10.*abs(co2_g_all_diffs[ppi,0,*,0,ki,z])) ; max_co2[ki] = max(10.*abs(co2_g_all_diffs[ppi,0,*,0,ki,z])) ; min_tot[ki] = min(diff_tot_source[ppi,0,*,0,ki,z]) ; max_tot[ki] = max(diff_tot_source[ppi,0,*,0,ki,z]) ;endfor ;poly = polygon([Kdiff,reverse(Kdiff)],[min_tot,reverse(max_tot)],/data,/fill_background,fill_color='medium orchid',fill_transparency=30,color='medium orchid',thick=2) ;poly = polygon([Kdiff,reverse(Kdiff)],[min_sspb,reverse(max_sspb)],/data,/fill_background,fill_color='aqua',fill_transparency=10,color='aqua',thick=2,pattern_orientation=45,pattern_spacing=2) ;poly = polygon([Kdiff,reverse(Kdiff)],[min_co2,reverse(max_co2)],/data,/fill_background,fill_color='yellow',fill_transparency=60,color='yellow',thick=2,pattern_orientation=45, $ ; pattern_spacing=2) ;;for pp = tr_peak[z-1]-1, tr_peak[z-1] do for ki = ki_prof[z]-30, ng-1 do p = plot(replicate(ki,npeaks_source),reform(abs(sspb_source[pp,0,*,0,ki,z])),symbol='o',/sym_filled,sym_size=0.5, $ ;; color='sky blue',linestyle=6,/overplot) ;;for pp = tr_peak[z-1]-1, tr_peak[z-1] do for ki = ki_prof[z]-30, ng-1 do p = plot(replicate(ki,npeaks_source),10.*reform(abs(co2_g_all_diffs[pp,0,*,0,ki,z])),symbol='o',/sym_filled,sym_size=0.5, $ ;; color='lime green',linestyle=6,/overplot) ;;for pp = tr_peak[z-1]-1, tr_peak[z-1] do for ki = ki_prof[z]-30, ng-1 do p = plot(replicate(ki,npeaks_source),reform(diff_tot_source[pp,0,*,0,ki,z]),symbol='o',/sym_filled,sym_size=0.5, $ ;; color=colors3[pp-1],linestyle=6,/overplot) ;if z ge 7 then begin ; pp1 = tr_peak[z-1]-1 ; pp2 = tr_peak[z-1] ;endif else begin ; pp1 = 0 ; pp2 = npeaks_source-1 ;endelse ;p = plot(replicate(Kdiff[ki_prof[z]],2),[0,40],linestyle=2,thick=2,/overplot) ;p = plot(replicate(Kdiff[ki_prof[z-1]],2),[0,40],linestyle=2,thick=2,color='grey',/overplot) ;p = plot(replicate(Kdiff[ki_prof[z-1]-3],2),[0,40],linestyle=1,thick=2,color='grey',/overplot) ;;s = symbol(Kdiff[ki_prof[z]],diff_tot_source[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_prof[z],z],symbol='Star',sym_size=1.5,sym_color=colors3[tr_peak[z]-1],/data) ;;for pp = pp1, pp2 do p = plot(replicate(Kdiff[ki_prof[z]],npeaks_source),reform(diff_tot_source[pp,0,*,0,ki_prof[z],z]),symbol='o',sym_size=1,color=colors3[pp-1],linestyle=6,/overplot) ;for pp = pp1, pp2 do for ppn = 0, npeaks_source-1 do s = symbol(Kdiff[ki_prof[z]],diff_tot_source[pp,0,ppn,0,ki_prof[z],z],'o',sym_size=ppn*0.25+1.0,sym_thick=2,sym_color=colors3[pp],/data) ;t = text(0.3,0.84,'D',color='medium orchid',font_size=12,font_style=2,/norm) ;t = text(0.3,0.8,'aSSPB$_{\mu}$',color='turquoise',font_size=12,/norm) ;t = text(0.3,0.76,'10*a$\chi_{CO2diff}$',color='gold',font_size=12,/norm) ;t = text(Kdiff[ki_prof[z]]+2,12,'K$_{opt}$(380K)',color='black',font_size=11,/data) ;t = text(Kdiff[ki_prof[z-1]]+4,11,'K$_{opt}$(370K)',color='grey',font_size=11,/data) ;t = text(0.28,0.57,'Tropical Lat Peak',font_size=11,/norm) ;t = text(0.35,0.53,strcompress(string(fix(tr_peaks[ppi[0]])),/r)+'$\deg$N',color=colors3[ppi[0]],font_size=11,font_style=1,/norm) ;t = text(0.35,0.49,strcompress(string(fix(tr_peaks[ppi[1]])),/r)+'$\deg$N',color=colors3[ppi[1]],font_size=11,font_style=1,/norm) ;t = text(0.31,0.44,'NH Lat Peak',font_size=11,/norm) ;s = symbol(0.34,0.41,'o',sym_size=(npeaks_source-1)*0.25+1.0,sym_thick=2,sym_color='black',/norm) ;s = symbol(0.34,0.37,'o',sym_size=1.0,sym_thick=2,sym_color='black',/norm) ;t = text(0.36,0.4,strcompress(string(fix(nh_peaks[-1])),/r)+'$\deg$N',color='black',font_size=11,/norm) ;t = text(0.36,0.36,strcompress(string(fix(nh_peaks[0])),/r)+'$\deg$N',color='black',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_diffs_vs_K_'+th1+'K.png' ;p = plot(indgen(2),/nodata,xrange=[10,1e3],/xlog,yrange=[0,3],xtitle='Diffusion coefficient (K,m$^2$/s)',ytitle='Differences', $ ; title='Minimum Combined and CO$_2$ Differences vs. Diffusion Coefficient',dimensions=[750,500],margin=[0.09,0.09,0.13,0.09],font_size=11) ;diffs = 0. & cdiffs = 0. & diffs_z = fltarr(nzth) & cdiffs_z = diffs_z ;for z = 6, 18 do begin ; restore,dir+'Aircraft/Missions/seac4rs/co2_predef_'+shapes[shp_i_opt[z]]+'.sav' ; restore,dir+'Aircraft/Missions/seac4rs/bl_predef_'+shapes[shp_i_opt[z]]+'.sav' ; p = plot(replicate(Kdiff[ki_prof[z]],2),[0,3],linestyle=2,thick=1,color='grey',/overplot) ; if z ge 7 then begin ; pp1 = tr_peak[z-1]-1 ; pp2 = tr_peak[z-1] ; endif else begin ; pp1 = 0 ; pp2 = npeaks_source-1 ; endelse ; diffs_z[z] = min(diff_tot_source[*,0,*,0,ki_prof[z],z]) ; cdiffs_z[z] = 10.*abs(co2_g_all_diffs[tr_peak[z],0,nh_peak[z],0,ki_prof[z],z]) ; diffs += min(diff_tot_source[*,0,*,0,ki_prof[z],z]) ; cdiffs += 10.*abs(co2_g_all_diffs[tr_peak[z],0,nh_peak[z],0,ki_prof[z],z]) ; for pp = pp1, pp2 do for ppn = 0, npeaks_source-1 do s = symbol(Kdiff[ki_prof[z]],diff_tot_source[pp,0,ppn,0,ki_prof[z],z],'o',sym_size=ppn*0.25+1.0,sym_thick=2,sym_color=colors3[pp],/data) ; s = symbol(Kdiff[ki_prof[z]],10.*abs(co2_g_all_diffs[tr_peak[z],0,nh_peak[z],0,ki_prof[z],z]),'s',sym_thick=2,sym_size=1.5,sym_color='medium orchid',/data) ; th1 = strcompress(string(fix(theta_grid[z])),/r) ; print,th1,min(diff_tot_source[*,0,*,0,ki_prof[z],z]),10.*abs(co2_g_all_diffs[tr_peak[z],0,nh_peak[z],0,ki_prof[z],z]) ; if z mod 2 eq 0 then offset = 0.1 else offset = 0 ; t = text(Kdiff[ki_prof[z]],2.8-offset,th1+'K',font_size=10,alignment=0.5,color='black',/data) ;endfor ;print,diffs,cdiffs ;;save,diffs_z,cdiffs_z,diffs,cdiffs,filename=dir+'Aircraft/Missions/seac4rs/diffs_'+shapes[shp_i]+'.sav' ;t = text(0.94,0.84,'Tropical',font_size=11,alignment=0.5,/norm) ;t = text(0.94,0.8,'Lat Peak',font_size=11,alignment=0.5,/norm) ;for y = 1, npeaks_source-1 do t = text(0.94,0.76-(y-1.)*0.04,strcompress(string(fix(tr_peaks[y])),/r)+'$\deg$N',color=colors3[y],font_size=11,alignment=0.5,/norm) ;t = text(0.94,0.44,'NH',font_size=11,alignment=0.5,/norm) ;t = text(0.94,0.4,'Lat Peak',font_size=11,alignment=0.5,/norm) ;s = symbol(0.91,0.36,'o',sym_size=(npeaks_source-1)*0.25+1.0,sym_thick=2,sym_color='black',/norm) ;s = symbol(0.91,0.32,'o',sym_size=1.0,sym_thick=2,sym_color='black',/norm) ;t = text(0.96,0.35,strcompress(string(fix(nh_peaks[-1])),/r)+'$\deg$N',color='black',font_size=11,alignment=0.5,/norm) ;t = text(0.96,0.31,strcompress(string(fix(nh_peaks[0])),/r)+'$\deg$N',color='black',font_size=11,alignment=0.5,/norm) ;t = text(0.94,0.23,'10*a$\chi_{CO2diff}$',font_size=11,alignment=0.5,color='medium orchid',/norm) ;s = symbol(0.94,0.19,'s',sym_size=1.5,sym_thick=2,sym_color='medium orchid',/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_min_diffs_vs_K.png' ggi = [2,3,4,5,6,7] ;ggi = indgen(6) diffs_z_all = fltarr(nzth,6) & cdiffs_z_all = diffs_z_all for gg = 0, 5 do begin restore,filename=dir+'Aircraft/Missions/seac4rs/diffs_'+shapes[ggi[gg]]+'.sav' diffs_z_all[*,gg] = diffs_z cdiffs_z_all[*,gg] = cdiffs_z endfor ;p = plot(indgen(2),/nodata,xrange=[0,3],yrange=[350,480],xtitle='aSSPB$_{\mu}$',ytitle='Theta (K)',title='Minimum Tracer Differences for Partitioning Choices',dimensions=[450,500], $ ; margin=[0.15,0.08,0.03,0.07],font_size=11) ;for gg = 0, 5 do p = plot(diffs_z_all[6:18,gg],theta_grid[6:18],thick=3,color=colors3[gg],symbol='o',/sym_filled,sym_size=1.5,/overplot) ;t = text(0.75,0.82,'30 days',color=colors3[0],font_size=11,/norm) ;t = text(0.75,0.78,'50 days',color=colors3[1],font_size=11,/norm) ;t = text(0.75,0.74,'100 days',color=colors3[2],font_size=11,/norm) ;t = text(0.75,0.7,'150 days',color=colors3[3],font_size=11,/norm) ;t = text(0.75,0.66,'200 days',color=colors3[4],font_size=11,/norm) ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_SSPB_min_profiles.png' ;; ;p = plot(indgen(2),/nodata,xrange=[0,3],yrange=[350,480],xtitle='10*a$\chi_{CO2diff}$',ytitle='Theta (K)',title='Minimum CO$_2$ Differences for Partitioning Choices',dimensions=[450,500], $ ; margin=[0.15,0.08,0.03,0.07],font_size=11) ;for gg = 0, 5 do p = plot(cdiffs_z_all[6:18,gg],theta_grid[6:18],thick=3,color=colors3[gg],symbol='o',/sym_filled,sym_size=1.5,linestyle=2,/overplot) ;t = text(0.75,0.82,'30 days',color=colors3[0],font_size=11,/norm) ;t = text(0.75,0.78,'50 days',color=colors3[1],font_size=11,/norm) ;t = text(0.75,0.74,'100 days',color=colors3[2],font_size=11,/norm) ;t = text(0.75,0.7,'150 days',color=colors3[3],font_size=11,/norm) ;t = text(0.75,0.66,'200 days',color=colors3[4],font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_min_diff_profiles.png' ;p = plot(indgen(2),/nodata,dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([ki_prof[z]-30,ng-1],[0,0],linestyle=1,/overplot) ;for pp = tr_peak[z-1]-1, tr_peak[z-1] do for ki = ki_prof[z]-30, ki_prof[z-1]-3, 3 do p = plot(replicate(ki,npeaks_source),reform(sspb_source[pp,0,*,0,ki,z]),symbol='o',/sym_filled,sym_size=0.5, $ ; color=colors3[pp-1],linestyle=6,/overplot) ;s = symbol(ki_prof[z],sspb_source[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],ki_prof[z],z],'s',/sym_filled,sym_color='blue',/data) x1 = [0,0.2,0.3,0.2,0.3,8e-2,2,1e-2,0.5,2,1,11,170,50,1.3,17,30,110,50,380,14] x2 = [700,5e3,2e3,5e3,1e3,50,1.5e4,10,40,17,10,29,270,100,4.4,26,90,250,77,540,17.5] xlog = [0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] units = ['ppb',replicate('ppt',20)] aa = 8 ;for i = 3, 13 do begin ; p = plot(indgen(2),/nodata,xrange=[x1[i],x2[i]],xlog=xlog[i],yrange=[290,490],xtitle=units[i],ytitle='Theta (K)',title='SEAC$^4$RS '+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4), $ ; dimensions=[400,500],margin=[0.15,0.08,0.03,0.07],font_size=11) ; p = plot(dc8_was[ci_vars_dc8[i]],dc8_was['THETA'],color='sky blue',symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ; p = plot(er2_was[ci_vars_er2[i]],er2_was['THETA'],color='blue',symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ;; p = plot(er2_was[ci_vars_er2[i]],er2_was['THETA'],rgb_table=39,vert_colors=(er2_was['LATITUDE']-15.)*8.,symbol='o',sym_size=0.5,linestyle=6,/sym_filled,/overplot) ;; chk = where(er2_was['THETA'] le er2_was['TropopausePotTemp1_MTP']) ;; p = plot(er2_was[ci_vars_er2[i],chk],er2_was['THETA',chk],color='lime green',symbol='o',sym_size=0.5,linestyle=6,/sym_filled,/overplot) ;; p = plot(er2_was_theta_grid[1,5:-1,where(varnames_er2_was eq ci_vars_er2[i])],theta_grid[5:-1],color='pink',thick=3,symbol='o',/sym_filled,/overplot) ;; p = plot(er2_was_theta_grid[0,5:-1,where(varnames_er2_was eq ci_vars_er2[i])],theta_grid[5:-1],color='medium orchid',thick=3,symbol='o',/sym_filled,/overplot) ; p = plot(dc8_was_theta_grid[0,0:4,where(varnames_dc8_was eq ci_vars_dc8[i])],theta_grid[0:4],color='dodger blue',thick=4,symbol='o',/sym_filled,sym_size=1.5,/overplot) ;; p = plot(ci_theta_grid[0,6:-1,i],theta_grid[6:-1],color='purple',thick=3,symbol='o',/sym_filled,errorbar_capsize=0,/overplot) ; gd = where(finite(ci_theta_grid[0,6:-1,i]),ngd) ; p = plot(ci_theta_grid[0,gd+6,i],theta_grid[gd+6],rgb_table=3,vert_colors=theta_grid[gd+6]-240,thick=4,linestyle=2,/overplot) ; p = plot(ci_theta_grid_adj[gd+6,i],theta_grid[gd+6],color='lime green',thick=3,/overplot) ; p = plot(ci_theta_grid_ind_adj[gd+6,i],theta_grid[gd+6],rgb_table=3,vert_colors=theta_grid[gd+6]-240,thick=4,symbol='o',/sym_filled,sym_size=1.5,/overplot) ;; s = symbol(er2_was[ci_vars_er2[i],zti[aa]],er2_was['THETA',zti[aa]],'s',sym_color='lime green',/sym_filled,/data) ;; p = plot(replicate(bl2[0,i],2),[theta_grid[0],theta_grid[2]],color='red',thick=3,/overplot) ;; p = plot(bl_best_source[gd+6,i],replicate(theta_grid[1],ngd),rgb_table=3,vert_colors=theta_grid[gd+6]-240,linestyle=6,symbol='o',/sym_filled,/overplot) ; p = plot(bl_best_source[gd+6,i],replicate(theta_grid[1],ngd),rgb_table=3,vert_colors=theta_grid[gd+6]-240,linestyle=6,symbol='o',/sym_filled,sym_size=1.5,/overplot) ; t = text(0.21,0.35,'ER2',color='blue',font_size=10,/norm) ;; t = text(0.76,0.8,'Avg',color='purple',/norm) ; t = text(0.21,0.31,'DC8',color='sky blue',font_size=10,/norm) ;; t = text(0.76,0.76,'Avg',color='dodger blue',/norm) ;; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'_theta_profile.png' ;; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+ci_vars_er2[i]+'_theta_profile1.png' ;; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+ci_vars_er2[i]+'_theta_profile2.png' ;endfor x1 = [0,0.4,0.7,1,6,1,1e2,8e-2,2,2,1,11,170,50,1,17,30,110,50,380,14] x2 = [700,5e3,9e2,1.5e4,7e2,80,1.5e4,30,40,10,17,29,270,100,4.2,26,90,250,77,540,17.5] ;for i = 3, 6 do begin ; p = plot(indgen(2),/nodata,xrange=[x1[i],x2[i]],xlog=xlog[i],yrange=[290,405],xtitle='Mixing ratio',ytitle='Theta (K)',title=strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),dimensions=[500,550], $ ; margin=[0.13,0.08,0.03,0.07],font_size=11) ; p = plot(dc8_was[ci_vars_dc8[i]],dc8_was['THETA'],color='sky blue',symbol='o',sym_size=0.5,linestyle=6,/sym_filled,/overplot) ; p = plot(er2_was[ci_vars_er2[i]],er2_was['THETA'],color='blue',symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ; chk = where(bl_nh_scale[*,i] gt 1,nchk) ; chk2 = where(bl_nh_scale[*,i] lt 1,nchk2) ; if nchk gt 0 then p = plot(er2_was[ci_vars_er2[i],zti[chk]],theta_ind[chk],color='red',symbol='o',sym_size=1,linestyle=6,sym_thick=1,/overplot) ;; p = plot(er2_was[ci_vars_er2[i],zti[chk2[0:2]]],theta_ind[chk2[0:2]],color='medium orchid',symbol='s',sym_size=1,sym_thick=2,linestyle=6,/overplot) ; if nchk2 gt 0 then p = plot(er2_was[ci_vars_er2[i],zti[chk2]],theta_ind[chk2],color='lime green',symbol='o',sym_size=1,linestyle=6,sym_thick=1,/overplot) ; p = plot(ci_theta_grid_adj[6:10,i],theta_grid[6:10],rgb_table=3,vert_colors=theta_grid[6:10]-240,thick=3,symbol='o',/sym_filled,/overplot) ;; p = plot(ci_theta_grid_ind_adj[6:10,i],theta_grid[6:10],color='gold',thick=3,symbol='o',/sym_filled,/overplot) ; p = plot(bl_nh_best_source[6:10,i]/(1.0-tr_best_source_frac[6:10]),replicate(theta_grid[0],ngd),rgb_table=3,vert_colors=theta_grid[6:10]-240,linestyle=6,symbol='o',/sym_filled,/overplot) ; if nchk gt 0 then p = plot(bl_nh_ind_scaled[chk,i],replicate(theta_grid[1],nchk),color='red',linestyle=6,symbol='o',sym_thick=2,/overplot) ;; p = plot(bl_nh_ind_scaled[chk2[0:2],i],replicate(theta_grid[1],3),color='medium orchid',linestyle=6,symbol='s',sym_thick=2,/overplot) ; if nchk2 gt 0 then p = plot(bl_nh_ind_scaled[chk2,i],replicate(theta_grid[1],nchk2),color='lime green',linestyle=6,symbol='o',sym_thick=2,/overplot) ; t = text(0.85,0.55,'ER2',color='blue',font_size=10,/norm) ; t = text(0.85,0.48,'DC8',color='sky blue',font_size=10,/norm) ; t = text(0.68,0.82,'Polluted NH source',color='red',font_size=10,/norm) ;; t = text(0.68,0.78,'Clean NH source',color='lime green',font_size=10,/norm) ; t = text(0.16,0.17,'Adjusted NH BLs',color='black',font_size=10,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'_theta_profile_NH_BL_adjust.png' ;endfor colors2 = ['dark slate gray','saddle brown','chocolate','peru','red','orange red','orange','goldenrod','gold','dark green','chartreuse','green','lime','lime green','aqua', $ 'yellow green','deep sky blue','dodger blue','blue','dark orchid','dark magenta','purple','magenta','plum'] colors4 = ['grey','orange red','orange','red','dodger blue','blue','gold','goldenrod','deep sky blue','peru','medium orchid','purple','lime green','green','yellow green'] ;p = plot(indgen(2),/nodata,xrange=[0,12],yrange=[340,485],xtitle='Normalized Mixing ratio',ytitle='Theta (K)',title='Normalized Mixing Ratios to Mean Theta Profiles', $ ; dimensions=[500,500],margin=[0.13,0.09,0.04,0.08],font_size=11) ;p = plot([1,1],[330,505],linestyle=2,/overplot) ;poly = polygon([0,12,12,0],[340,340,415,415],/fill_background,/data,linestyle=6,fill_color='orange',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;for i = 3, 5, 2 do begin ; p = plot(mean_prof_rel[*,i],theta_ind,color=colors4[i],symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ; t = text(0.65,0.86-(i-1)*0.02,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=10,/norm) ;endfor ;for i = 6, 6, 2 do begin ; p = plot(mean_prof_rel[*,i],theta_ind[*],color=colors4[i],symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ; t = text(0.65,0.84-(i-1)*0.02,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=10,/norm) ;endfor ;;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_enhanced_tracers_theta_profile1.png' ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_enhanced_tracers_theta_profile2.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_enhanced_tracers_theta_profile.png' ;p = plot(indgen(2),/nodata,xrange=[0.7,2],yrange=[330,485],xtitle='Normalized Mixing ratio',ytitle='Theta (K)',title='SEAC$^4$RS Normalized Mixing Ratios to Mean Theta Profile', $ ; dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([1,1],[330,505],linestyle=2,/overplot) ;for i = 9, 18 do begin ; p = plot(mean_prof_rel[*,i],theta_ind,color=colors2[i],symbol='o',sym_size=0.75,linestyle=6,/sym_filled,/overplot) ; t = text(0.7,0.84-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors2[i],/norm) ;endfor ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_enhanced_tracers_theta_profile3.png' p = plot(indgen(2),/nodata,xrange=[-0.02,0.75],yrange=[345,485],xtitle='Fraction',ytitle='Theta (K)',title='NH Extratropical BL Source Fraction',dimensions=[400,500], $ margin=[0.16,0.08,0.03,0.07],font_size=12) p = plot(1.-tr_source_frac,theta_ind,color='sky blue',symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) p = plot(1.-tr_best_source_frac[6:18],theta_grid[6:18],color='blue',symbol='s',thick=3,sym_size=1,/sym_filled,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_NH_BL_frac.png' i = 1 ;p = plot(indgen(2),/nodata,xrange=[0,11],yrange=[340,485],xtitle='Ratio',ytitle='Theta (K)',title='BL Fraction Ratios to Mean Theta Profile',dimensions=[500,500],font_size=11, $ ; margin=[0.13,0.09,0.04,0.08]) ;p = plot([1,1],[340,505],linestyle=2,/overplot) ;p = plot(tf_ind_init_mean_rel[*,i],theta_ind,color=colors4[i],symbol='o',sym_size=1.5,linestyle=6,/overplot) ;p = plot(tf_ind_best_mean_rel[*,i],theta_ind,color=colors4[i],symbol='s',sym_size=1,linestyle=6,/sym_filled,/overplot) ;t = text(0.7,0.84-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=10,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_ratio_profile_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.png' ; ;p = plot(indgen(2),/nodata,xrange=[5e-3,4.5],/xlog,yrange=[340,485],xtitle='$\mu$',ytitle='Theta (K)',title='BL Fractions',dimensions=[500,500],margin=[0.13,0.09,0.04,0.08],font_size=11) ;p = plot([1,1],[340,505],linestyle=2,/overplot) ;p = plot(tf_ind_init[*,i],theta_ind,color=colors4[i],symbol='o',sym_size=1.5,linestyle=6,/overplot) ;p = plot(tf_best_ind[*,i],theta_ind,color=colors4[i],symbol='s',sym_size=1,linestyle=6,/sym_filled,/overplot) ;t = text(0.78,0.84-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=10,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_ratio_ind_profile_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.png' ; ;for i = 3, 3 do begin ; p = plot(indgen(2),/nodata,xrange=[0.7,100],/xlog,yrange=[345,420],xtitle='Scale Factor',ytitle='Theta (K)',title='NH Boundary Layer Scaling',dimensions=[450,550],margin=[0.14,0.08,0.03,0.07], $ ; font_size=11) ; p = plot([1,1],[330,505],linestyle=2,/overplot) ; gd = where(bl_nh_scale[*,i] ne 1) ; p = plot(bl_nh_scale[gd,i],theta_ind[gd],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ; t = text(0.5,0.82,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=11,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_nh_scale_profile_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.png' ; print,bl_nh_scale[gd,i];,reform(scale_tf_tau[i,*]),reform(scale_tf[i,*]) ;endfor ;p = plot(indgen(2),/nodata,xrange=[0.9,100],/xlog,yrange=[0.9,100],/ylog,xtitle='Propane Scale Factor',ytitle='Tracer Scale Factor',title='NH BL Scaling Correlations',dimensions=[700,500], $ ; margin=0.09,font_size=11) ;p = plot([0.3,1e2],[0.3,1e2],linestyle=2,/overplot) ;p = plot([1,1],[1,1e2],linestyle=2,/overplot) ;p = plot([1,1e2],[1,1],linestyle=2,/overplot) ;for i = 1, 2 do begin ; gd = where(bl_nh_scale[*,i] ne 1) ; p = plot(bl_nh_scale[gd,3],bl_nh_scale[gd,i],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;endfor ;gd = where(bl_nh_scale[*,3] ne 1) ;p = plot(bl_nh_scale[gd,3],bl_nh_scale[gd,3],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/overplot) ;for i = 4, 9 do begin ; gd = where(bl_nh_scale[*,i] ne 1) ; p = plot(bl_nh_scale[gd,3],bl_nh_scale[gd,i],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;endfor ;for i = 1, 9 do t = text(0.2,0.85-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_relative_to_propane.png' ;p = plot(indgen(2),/nodata,xrange=[0.6,70],/xlog,yrange=[0.6,50],/ylog,xtitle='Ethyne Scale Factor',ytitle='Tracer Scale Factor',title='NH BL Scaling Correlations',dimensions=[700,500], $ ; margin=0.09,font_size=11) ;p = plot([0.3,60],[0.3,60],linestyle=2,/overplot) ;p = plot([1,1],[1,60],linestyle=2,/overplot) ;p = plot([1,60],[1,1],linestyle=2,/overplot) ;for i = 1, 3 do begin ; gd = where(bl_nh_scale[*,i] ne 1) ; p = plot(bl_nh_scale[gd,4],bl_nh_scale[gd,i],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;endfor ;for i = 4, 9 do begin ; gd = where(bl_nh_scale[*,i] ne 1) ; p = plot(bl_nh_scale[gd,4],bl_nh_scale[gd,i],color=colors4[i],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;endfor ;for i = 1, 9 do t = text(0.21,0.84-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],/norm) z = 8 ;p = plot(indgen(2),/nodata,xrange=[0,1.1],yrange=[0,1]) ;p = plot([1,1],[0,6],linestyle=2,/overplot) ;;p = plot(H_i,hist_all_mean_rel[*,i],thick=1,color=colors2[i],/overplot) ;for i = 1, nciv-1 do begin ;; hsm = hist_tf_ind_init_prof[*,z,i] ;; hsm /= max(hsm,mi) ;; p = plot(H_i,hsm,thick=2,color=colors2[i],/overplot) ; hsm = hist_tf_ind_adj_prof[*,z,i] ; hsm /= max(hsm,mi) ; p = plot(H_i,hsm,thick=2,color=colors2[i],/overplot) ;; hsm = smooth(hist_tf_ind_init_mean_rel[*,i],11,/edge_truncate) ;; hsm /= max(hsm,mi) ;; p = plot(H_i,hsm,thick=2,color=colors2[i],linestyle=1,/overplot) ;; for z = 6, 6 do begin ;; hsm = hist_tf_ind_init_mean_rel_prof[*,z,i] ;; hsm /= max(hsm,mi) ;; p = plot(H_i,hsm,thick=3,color=colors2[i],/overplot) ;; endfor ;; chk = where(hsm ge 0.5) ;; half = 0.1*(H_i[chk[-1]] - H_i[chk[0]]) ;; gauss = exp(-(H_i-H_i[mi])^2. / (2.*half)) ;; gauss /= max(gauss) ;; p = plot(H_i,gauss,thick=3,color=colors2[i],linestyle=2,/overplot) ;; hsm = smooth(hist_tf_ind_adj_mean_rel[*,i],11,/edge_truncate) ;; hsm /= max(hsm,mi) ;; p = plot(H_i,hsm,thick=3,color=colors2[i],linestyle=1,/overplot) ;; hsm = smooth(hist_all_mean_rel[*,i],11,/edge_truncate) ;; hsm /= max(hsm,mi) ;; p = plot(H_i,hsm,thick=3,color=colors2[i],linestyle=1,/overplot) ;endfor ;p = plot(indgen(2),/nodata,xrange=[0.4,1.6]) ;p = plot([1,1],[0,1],linestyle=2,/overplot) ;for i = 8, 20 do begin ; hsm = hist_all_mean_rel[*,i] ; p = plot(H_i,hsm/max(hsm),thick=3,color=colors2[i],/overplot) ;endfor i = 12 ;p = plot(indgen(2),/nodata,yrange=[0.3,2]) ;gd = where(mean_prof_rel[*,i] ne 1) ;p = plot(lat_ind[gd],mean_prof_rel[gd,i],symbol='o',color='blue',/sym_filled,linestyle=6,/overplot) ;for z = 6, 6 do begin ; zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;endfor ;p = plot(indgen(2),/nodata,yrange=[0,3.5],xrange=[0,11],xtitle='Normalized Propane',ytitle='Normalized Mixing Ratio',title='SEAC4RS Normalized Mixing Ratios to Mean Theta Profile') ;p = plot([1,1],[0,12],linestyle=2,/overplot) ;p = plot([0,12],[1,1],linestyle=2,/overplot) ;for i = 5, 8, 3 do begin ; p = plot(mean_prof_rel[*,3],mean_prof_rel[*,i],color=colors[i-1],symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ; chk = where(finite(mean_prof_rel[*,i]) and mean_prof_rel[*,i] gt 0 and finite(mean_prof_rel[*,3]) and mean_prof_rel[*,3] gt 0) ;; fit = linfit(mean_prof_rel[chk,i],mean_prof_rel[chk,3]) ;; print,fit,fit[0] + fit[1]*[0,7] ;; p = plot([0,7],fit[0] + fit[1]*[0,7],color=colors[i],/overplot) ;endfor ;p = plot(indgen(2),/nodata,xrange=[0,ng],yrange=[-5,4],xtitle='Ki',ytitle='ppmv',title='SEAC4RS CO2 Age Derived Diffs') ;p = plot([0,ng],[0,0],linestyle=1,/overplot) ;for z = 10, 10 do begin ; for y = 7, ny_ct-4 do p = plot(indgen(ng),reform(co2_g[y,*,z])-(er2_was_theta_grid[1,z,where(varnames_er2_was eq 'CO2_HUPCRS')])[0],thick=2,color=colors[y-6],/overplot) ; p = plot(replicate(diff_min_i[z],ny_ct),reform(co2_g[*,diff_min_i[z],z])-(er2_was_theta_grid[1,z,where(varnames_er2_was eq 'CO2_HUPCRS')])[0],symbol='o',/sym_filled,color='medium orchid', $ ; linestyle=6,/overplot) ;endfor months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec','Jan'] ;for i = 6, 6 do begin ; tti = where(co2_time_ext ge 2012.5 and co2_time_ext le 2014.5) ; print,flask_data_ext[tti,2:3,i] ; mn = 0.9*min(flask_data_ext[tti,1:3,i]) ; mm = 1.1*max(flask_data_ext[tti,1:3,i]) ; p = plot(indgen(2),/nodata,xrange=[2013,2014],xtickname=months,xminor=0,xmajor=13,xtitle='2013',yrange=[mn,mm],ytitle='Mixing Ratio',dimensions=[700,350],margin=0.1,font_size=11, $ ; title='GML NH Extratropical Surface '+strmid(ci_vars_er2[i],0,strlen(ci_vars_er2[i])-4)+' Time Series') ; poly = polygon([min(er2_yrfrac),max(er2_yrfrac),max(er2_yrfrac),min(er2_yrfrac)],[0,0,mm,mm],/fill_background,/data,linestyle=6,fill_color='red',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ; p = plot([er2_yrfrac_mean,er2_yrfrac_mean],[0,mm],linestyle=2,thick=3,color='red',/overplot) ;; for y = 40, 90, 5 do p = plot(co2_time_ext,surface_data_interp[*,y,i],thick=3,color='medium orchid',/overplot) ;; for pp = 0, 3 do p = plot(co2_time_ext,surf[*,pp,0,1,i],thick=3,color='red',/overplot) ; for y = 1, 3 do p = plot(co2_time_ext,flask_data_ext[*,y,i],thick=3,color=colors[y],/overplot) ;; for y = 0, 5 do p = plot(co2_time,flask_data[*,y,i],thick=2,color=colors[y],/overplot) ;; t = text(0.84,0.23,'RPB',color=colors[1],/norm) ;; t = text(0.84,0.26,'THD',color=colors[2],/norm) ;; t = text(0.84,0.29,'MHD',color=colors[3],/norm) ;; t = text(0.84,0.29,'LEF',color=colors[3],/norm) ;; t = text(0.84,0.26,'NWR',color=colors[2],/norm) ;; t = text(0.84,0.23,'KEY',color=colors[1],/norm) ; t = text(0.2,0.28,'45$\deg$-60$\deg$N',color=colors[3],font_size=11,/norm) ; t = text(0.2,0.24,'20$\deg$-45$\deg$N',color=colors[2],font_size=11,/norm) ; t = text(0.2,0.2,'0$\deg$-20$\deg$N',color=colors[1],font_size=11,/norm) ;; t = text(0.84,0.23,'0-20N',color=colors[1],/norm) ;; t = text(0.84,0.23,'Tropics',color=colors[4],/norm) ;; t = text(0.84,0.2,'NH',color=colors[5],/norm) ; t = text(0.7,0.8,'SEAC$^4$RS',color='red',font_size=10,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+strmid(ci_vars_er2[i],0,strlen(ci_vars_er2[i])-4)+'_NH_time_series.png' ;endfor ;for i = 6, 6 do begin ; mn = min(flask_data_ext[*,0:1,i]) ; mm = max(flask_data_ext[*,0:1,i]) ; p = plot(indgen(2),/nodata,xrange=[2003,2014],xtitle='Year',yrange=[mn,mm],ytitle='Mixing Ratio',dimensions=[700,350],margin=0.1,font_size=11, $ ; title='GML Tropical Surface '+strmid(ci_vars_er2[i],0,strlen(ci_vars_er2[i])-4)+' Time Series') ; poly = polygon([min(er2_yrfrac),max(er2_yrfrac),max(er2_yrfrac),min(er2_yrfrac)],[0,0,mm,mm],/fill_background,/data,linestyle=6,fill_color='red',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ; p = plot([er2_yrfrac_mean,er2_yrfrac_mean],[0,mm],linestyle=2,thick=2,color='red',/overplot) ;; for pp = 0, 3 do p = plot(co2_time_ext,surf[*,pp,0,1,i],thick=3,color='red',/overplot) ;; for pp = 0, npeaks_source-1 do p = plot(co2_time_ext,surf[*,pp,0,1,i],thick=2,color='red',/overplot) ;; for pp = 0, npeaks_source-1 do p = plot(co2_time_ext,surf[*,pp,1,1,i],thick=2,color='orange',/overplot) ;; for pp = 0, npeaks_source-1 do p = plot(co2_time_ext,surf[*,pp,3,1,i],thick=2,color='lime green',/overplot) ; for y = 0, 1 do p = plot(co2_time_ext,flask_data_ext[*,y,i],thick=2,color=colors[y],/overplot) ;; for y = 0, 30, 5 do p = plot(co2_time_ext,surface_data_interp[*,y,i],thick=2,color='medium orchid',/overplot) ; ; for y = 0, 5 do p = plot(co2_time,flask_data[*,y,i],thick=2,color=colors[y],/overplot) ; t = text(0.2,0.82,'0$\deg$-20$\deg$N',color=colors[1],font_size=11,/norm) ; t = text(0.2,0.78,'20$\deg$S-0$\deg$',color=colors[0],font_size=11,/norm) ; ; t = text(0.84,0.23,'0-20N',color=colors[1],/norm) ;; t = text(0.2,0.82,'RPB',color=colors[1],/norm) ;; t = text(0.2,0.78,'SMO',color=colors[0],/norm) ; ; t = text(0.84,0.23,'Tropics',color=colors[4],/norm) ; ; t = text(0.84,0.2,'NH',color=colors[5],/norm) ; t = text(0.78,0.8,'SEAC$^4$RS',color='red',font_size=10,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_'+strmid(ci_vars_er2[i],0,strlen(ci_vars_er2[i])-4)+'_tropical_time_series.png' ;endfor ; ;p = plot(indgen(2),/nodata,xrange=[2013,2014],yrange=[0.9,1.2],xtitle='Year',ytitle='Normalized Mixing Ratio',title='GML NH Surface Time Series') ;p = plot([er2_yrfrac_mean,er2_yrfrac_mean],[0,10],linestyle=2,/overplot) ;p = plot([2013,2014],[1,1],linestyle=2,/overplot) ;for i = 8, 8 do p = plot(co2_time_ext,flask_data_ext[*,5,i]/mean(flask_data_ext[-63:-62,5,i]),thick=2,color=colors[i-1],/overplot) ; ;p = plot(indgen(2),/nodata,xrange=[1983,2016],yrange=[0.1,1.4],xtitle='Year',ytitle='Normalized Mixing Ratio',title='GML Tropical Surface Time Series') ;p = plot([er2_yrfrac_mean,er2_yrfrac_mean],[0,3],linestyle=2,/overplot) ;p = plot([1983,2014],[1,1],linestyle=2,/overplot) ;for i = 9, 12 do p = plot(co2_time_ext,flask_data_ext[*,4,i]/mean(flask_data_ext[-68:-56,4,i]),thick=2,color=colors[i-1],/overplot) cgloadct,43,rgb_table=rgb_43 ;p = plot(indgen(2),/nodata,xrange=[0,1.02],xtitle='Boundary Layer Fraction ($\mu$)',yrange=[350,500],ytitle='Theta (K)',title='SEAC$^4$RS BL Fractions',dimensions=[700,500], $ ; margin=[0.08,0.09,0.18,0.08]) ;p = plot([1,1],[345,485],linestyle=2,/overplot) ;for i = 1, 20 do begin ;; p = plot(tf_best_source[*,i],theta_grid,color=colors2[i],thick=2,symbol='o',/sym_filled,/overplot) ; p = plot(tf_profile_ind[*,i],theta_grid,color=colors2[i],thick=2,symbol='o',/sym_filled,/overplot) ; t = text(0.83,0.75-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors2[i],/norm) ;endfor ;;p = plot(trop_fracs_ind_gradj[-1,aa,*],replicate(interpolate(theta_grid,zi[aa]),nciv),color='pink',thick=2,symbol='o',/sym_filled,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_profiles_theta.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_profiles_source_theta.png' ;for i = 6, 6, 3 do begin ; p = plot(indgen(2),/nodata,xtitle='$\mu$',xrange=[5e-2,5],/xlog,yrange=[340,420],ytitle='Theta (K)',title='Boundary Layer Fractions',dimensions=[450,550],margin=[0.14,0.08,0.03,0.07],font_size=11) ; p = plot([1,1],[345,485],linestyle=2,/overplot) ; p = plot(tf_best_source[*,i],theta_grid,color=colors4[i],thick=4,symbol='o',/sym_filled,/overplot) ; ;p = plot(tf_profile_ind[*,i],theta_grid,color=colors2[i],thick=3,symbol='o',/sym_filled,/overplot) ; p = plot(tf_ind_init[*,i],theta_ind,color=colors4[i],symbol='o',sym_size=1.5,linestyle=6,/overplot) ; ;p = plot(tf_ind_adj[*,i],theta_ind,color=colors2[i],symbol='o',sym_size=0.75,/sym_filled,linestyle=6,/overplot) ; p = plot(tf_best_ind[*,i],theta_ind,color=colors4[i],symbol='s',sym_size=1,/sym_filled,linestyle=6,transparency=20,/overplot) ; ;p = plot(tf_ind_adj_mean[*,i],theta_grid,color='medium orchid',thick=3,symbol='o',/sym_filled,/overplot) ; t = text(0.81,0.75,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=11,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_ind_profile_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.png' ;endfor ;p = plot(indgen(2),/nodata,xtitle='$\tau$ (Days)',xrange=[1,1e5],/xlog,yrange=[345,480],ytitle='Theta (K)',title='Trace Gas Lifetimes',dimensions=[700,500],margin=0.09,font_size=11) ;for i = 1, 20 do begin ; p = plot(tau_profile_source[i,6:18],theta_grid[6:18],color=colors2[i],thick=2,symbol='o',/sym_filled,/overplot) ; p = plot(tau_profile_source_init[i,6:18],theta_grid[6:18],color=colors2[i],thick=2,symbol='o',linestyle=2,/overplot) ; p = plot(tau_profile_source_inter[i,6:18],theta_grid[6:18],color=colors2[i],thick=2,symbol='o',linestyle=1,/overplot) ;; p = plot(tau_profile_ind[i,6:20],theta_grid[6:20],color=colors2[i],thick=2,symbol='o',/sym_filled,/overplot) ; t = text(0.13,0.85-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors2[i],/norm) ;endfor ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_lifetime_profiles_source_theta.png' ;for i = 6, 6, 3 do begin ; p = plot(indgen(2),/nodata,xtitle='$\tau$ (Days)',xrange=[10,4e2],/xlog,yrange=[345,420],ytitle='Theta (K)',title='Path Integrated Lifetimes',dimensions=[450,550],margin=[0.14,0.09,0.04,0.08], $ ; font_size=11) ; ;p = plot(tau_profile_source_init[i,6:18],theta_grid[6:18],color=colors2[i],thick=2,symbol='o',linestyle=2,/overplot) ; p = plot(tau_profile_source[i,6:18],theta_grid[6:18],color=colors4[i],thick=4,symbol='o',sym_size=2,sym_thick=2,/overplot) ; p = plot(0.25*tau_profile_source[i,6:18],theta_grid[6:18],color='grey',thick=2,linestyle=1,/overplot) ; p = plot(4.*tau_profile_source[i,6:18],theta_grid[6:18],color='grey',thick=2,linestyle=1,/overplot) ; ;p = plot(tau_profile_ind[i,6:20],theta_grid[6:20],color='blue',thick=2,symbol='o',/sym_filled,/overplot) ; ;p = plot(tau_ind_init[i,*],theta_ind,color='sky blue',symbol='o',/sym_filled,sym_size=0.5,linestyle=6,/overplot) ; gd = where(finite(tf_best_ind[*,i])) ; p = plot(tau_ind[i,gd],theta_ind[gd],color=colors4[i],symbol='s',/sym_filled,sym_size=1,linestyle=6,/overplot) ; t = text(0.78,0.67,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],font_size=11,/norm) ; p.save,dir+'Plots/SEAC4RS/SEAC4RS_tau_ind_profile_'+strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+'.png' ;endfor ;p = plot(indgen(2),/nodata,xtitle='ppmv',xrange=[-0.2,0.2],yrange=[350,505],ytitle='Theta (K)',title='CO$_2$ Differences',dimensions=[500,500],margin=[0.13,0.09,0.04,0.08],font_size=11) ;p = plot([0,0],[345,505],linestyle=2,/overplot) ;p = plot(co2_best_source-er2_was_theta_grid[1,*,where(varnames_er2_was eq 'CO2_HUPCRS')],theta_grid,color='purple',symbol='o',/sym_filled,linestyle=6,/overplot) ;p = plot(co2_best_ind-co2_ind,theta_ind,color='sky blue',symbol='o',/sym_filled,sym_size=0.5,linestyle=6,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_CO2_diff_profile.png' ;p = plot(indgen(2),/nodata,xtitle='Years',yrange=[345,485],xrange=[0,4],ytitle='Theta (K)',title='Mean Ages',dimensions=[400,500],margin=[0.16,0.08,0.03,0.07],font_size=12) ;;p = plot(mean_age_theta/365.,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,/overplot) ;p = plot(mean_age_ind/365.,theta_ind,color='sky blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(mean_age_source/365.,theta_grid,color='blue',thick=3,symbol='s',sym_size=1,/sym_filled,/overplot) ;;p = plot(mean_age_ind_tr/365.,theta_ind,color='orange',symbol='o',/sym_filled,sym_size=0.5,linestyle=6,/overplot) ;;p = plot(mean_age_source_tr/365.,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,linestyle=2,/overplot) ;;p = plot(mean_age_source_nh/365.,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,linestyle=1,/overplot) ;t = text(0.45,0.26,'Mean profile',color='blue',font_size=11,/norm) ;t = text(0.45,0.22,'Individual measurements',color='sky blue',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_mean_age_profile_theta.png' p = plot(indgen(2),/nodata,xtitle='Days',yrange=[345,485],xrange=[0.8,500],/xlog,ytitle='Theta (K)',title='Modal Ages',dimensions=[400,500],margin=[0.16,0.08,0.03,0.07],font_size=12) ;p = plot(modal_age_ind_tr,theta_ind,color='orange',symbol='o',/sym_filled,sym_size=0.5,linestyle=6,/overplot) p = plot(modal_age_ind,theta_ind,color='sky blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) p = plot(modal_age_source,theta_grid,color='blue',thick=3,symbol='s',/sym_filled,sym_size=1,/overplot) ;chk = where(modal_age_source_tr eq modal_age_source) ;p = plot(modal_age_source[chk],theta_grid[chk],color='lime green',symbol='o',sym_size=1.5,sym_thick=2,linestyle=6,/overplot) ;chk = where(modal_age_source_nh eq modal_age_source) ;p = plot(modal_age_source[chk],theta_grid[chk],color='medium orchid',symbol='o',sym_size=1.5,sym_thick=2,linestyle=6,/overplot) ;chk = where(modal_age_ind_tr eq modal_age_ind) ;p = plot(modal_age_ind[chk],theta_ind[chk],color='lime green',symbol='o',sym_size=1,sym_thick=2,linestyle=6,/overplot) ;chk = where(modal_age_ind_nh eq modal_age_ind) ;p = plot(modal_age_ind[chk],theta_ind[chk],color='medium orchid',symbol='o',sym_size=1,sym_thick=2,linestyle=6,/overplot) ;p = plot(modal_age_source_tr,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,linestyle=2,/overplot) ;p = plot(modal_age_source_nh,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,linestyle=1,/overplot) ;t = text(0.6,0.24,'Tropical peak',color='lime green',font_size=11,/norm) ;t = text(0.6,0.2,'NH peak',color='medium orchid',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_modal_age_profile_theta.png' ;p = plot(indgen(2),/nodata,xtitle='Years',yrange=[350,480],ytitle='Theta (K)',title='SEAC$^4$RS Mean and Modal Ages',dimensions=[500,600],margin=[0.13,0.15,0.04,0.08],font_size=11) ;p = plot(mean_age_source/365.,theta_grid,color='blue',thick=3,symbol='o',/sym_filled,/overplot) ;p = plot(indgen(2),/nodata,xtitle='Latitude',xrange=[-10,60],yrange=[345,485],ytitle='Theta (K)',title='BL Source Peak Latitudes',dimensions=[450,550],margin=[0.14,0.08,0.03,0.07],font_size=11) ;p = plot(tr_peaks[tr_peak[6:18]],theta_grid[6:18],color='blue',thick=3,/overplot) ;for z = 6, 18 do s = symbol(tr_peaks[tr_peak[z]],theta_grid[z],sym_color='blue','o',/sym_filled,sym_size=tr_best_source_frac[z]*2.5,/data) ;p = plot(nh_peaks[nh_peak[6:12]],theta_grid[6:12],color='blue',thick=3,/overplot) ;for z = 6, 12 do s = symbol(nh_peaks[nh_peak[z]],theta_grid[z],sym_color='blue','o',/sym_filled,sym_size=(1.-tr_best_source_frac[z])*2.5,/data) ;;p = plot(tr_source_profile_ind[6:18],theta_grid[6:18],color='green',thick=3,/overplot) ;;p = plot(nh_source_profile_ind[6:12],theta_grid[6:12],color='green',thick=3,/overplot) ;gd = where(tr_source_frac gt 0.2 and co2_ind gt 0,ngd) ;for aa = 0, ngd-1 do s = symbol(tr_peaks[tr_peak_ind[gd[aa]]],theta_ind[gd[aa]],sym_color='lime green','o',/sym_filled,sym_size=tr_source_frac[gd[aa]]*1.5,/data) ;gd = where(tr_source_frac lt 0.8 and co2_ind gt 0,ngd) ;for aa = 0, ngd-1 do s = symbol(nh_peaks[nh_peak_ind[gd[aa]]],theta_ind[gd[aa]],sym_color='medium orchid','o',/sym_filled,sym_size=(1.-tr_source_frac[gd[aa]])*1.5,/data) ;;p = plot(nh_peaks[nh_peak_ind[gd]],theta_ind[gd],color='medium orchid',symbol='o',/sym_filled,linestyle=6,sym_size=0.5,/overplot) ;t = text(0.6,0.82,'Tropical',color='lime green',font_size=11,/norm) ;t = text(0.6,0.78,'NH Extratropical',color='medium orchid',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_source_peak_lat_profiles.png' ;p = plot(indgen(2),/nodata,xtitle='Tracer Lifetime',xrange=[1,1e5],/xlog,yrange=[340,480],ytitle='Theta (K)',title='SEAC4RS Tracer Lifetimes') ;p = plot(ci_lifetimes_theta_profile[-1,14,*],theta_grid,color='purple',thick=2,symbol='o',/sym_filled,/overplot) ;p = plot(ci_lifetimes_ind[-1,*,14],theta_ind,color='pink',symbol='o',/sym_filled,linestyle=6,/overplot) ;p = plot(indgen(2),/nodata,xtitle='Modal Age (Days)',yrange=[0,500],xrange=[0.8,50],/xlog,ytitle='Ozone (ppbv)',title='Modal Ages vs Ozone',dimensions=[450,550],margin=[0.15,0.08,0.03,0.07], $ ; font_size=12) ;p = plot(modal_age_ind,o3_ind,color='sky blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(modal_age_source[6:12],ci_theta_grid[0,6:12,0],color='blue',symbol='o',/sym_filled,sym_size=1.5,thick=3,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_modal_age_vs_ozone.png' ;p = plot(indgen(2),/nodata,xtitle='Modal Age (Days)',yrange=[0,500],xrange=[0,50],ytitle='Ozone (ppbv)',title='Modal Ages vs Ozone',dimensions=[450,550],margin=[0.15,0.08,0.03,0.07], $ ; font_size=12) ;p = plot(modal_age_ind,o3_ind,color='sky blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(modal_age_source[6:12],ci_theta_grid[0,6:12,0],color='blue',symbol='o',/sym_filled,sym_size=1.5,thick=3,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_modal_age_vs_ozone_lin.png' color_mult = 175 ;p = plot(indgen(2),/nodata,ytitle='Theta (K)',xrange=[0,500],yrange=[345,405],xtitle='O3 (ppbv)',title='Ozone',dimensions=[500,500],margin=[0.16,0.08,0.03,0.07],font_size=12) ;mc = alog10(modal_age_ind)*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(o3_ind,theta_ind,rgb_table=39,vert_colors=mc,symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) i = 9 ;p = plot(indgen(2),/nodata,xrange=[0,55],yrange=[290,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='SEAC$^4$RS Sampling',margin=0.09,font_size=10.5, $ ; dimensions=[700,500]) ;p = plot(er2_was['LATITUDE'],er2_was['TropopausePotTemp1_MTP'],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) ;;p = plot(tpzm['lat','value',215:-1],tpzm_theta[215:-1,404],color='dark orange',thick=3,/overplot) ;poly = polygon([mer_lat[215:-1],reverse(mer_lat[215:-1])],[tpm_theta_avg[215:-1]-tpm_theta_sd_avg[215:-1],reverse(tpm_theta_avg[215:-1])+reverse(tpm_theta_sd_avg[215:-1])],/fill_background,/data, $ ; linestyle=6,fill_color='orange',transparency=60) ;p = plot(mer_lat[215:-1],tpm_theta_avg[215:-1],color='dark orange',thick=3,/overplot) ;;p = plot(mer_lat,tpm_theta2,color='gold',thick=3,/overplot) ;;p = plot(dc8_was['LATITUDE'],dc8_was['THETA'],color='sky blue',symbol='o',sym_size=0.5,linestyle=6,/sym_filled,/overplot) ;p = plot(er2_was['LATITUDE'],er2_was['THETA'],color='blue',symbol='o',sym_size=0.5,linestyle=6,/sym_filled,/overplot) ;;p = plot([20,24.5,36.5,46],[315,295,295,295],symbol='s',/sym_filled,color='red',linestyle=6,/overplot) ;;s = symbol(lat_ind[14],theta_ind[14],'s',sym_color='purple',/sym_filled,/data) ;;s = symbol(lat_ind[14],tpause_theta_ind[14],'s',sym_color='red',/sym_filled,/data) ;t = text(0.2,0.65,'WAS ER-2',color='blue',font_size=10,/norm) ;t = text(0.2,0.37,'Tropopause',color='orange',font_size=10,/norm) ;;t = text(0.2,0.22,'WAS DC-8',color='sky blue',font_size=10,/norm) ;;t = text(0.2,0.14,'NOAA sites',color='red',font_size=10,/norm) ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_sampling_theta.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_sampling_theta2.png' color_mult = 80. ;p = plot(indgen(2),/nodata,xrange=[15,53],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='SEAC$^4$RS Mean Ages',margin=0.105,font_size=11, $ ; dimensions=[650,400]) ;ii = where(finite(mean_age_ind)) ;;p = plot(lat_ind[ii],tpause_theta_ind[ii],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) ;poly = polygon([mer_lat[215:-1],reverse(mer_lat[215:-1])],[tpm_theta_avg[215:-1]-tpm_theta_sd_avg[215:-1],reverse(tpm_theta_avg[215:-1])+reverse(tpm_theta_sd_avg[215:-1])],/fill_background,/data, $ ; linestyle=6,fill_color='orange',transparency=70) ;p = plot(mer_lat[215:-1],tpm_theta_avg[215:-1],color='dark orange',thick=3,/overplot) ;mc = mean_age_ind[ii]/365.*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(lat_ind[ii],theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;tickvals = findgen(13)*0.25 ;tickname = strmid(strcompress(string(tickvals),/r),0,4) ;cb = COLORBAR(POSITION=[0.9015,0.23,0.9165,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='Years') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_mean_age_ind.png' ; color_mult = 100. ;p = plot(indgen(2),/nodata,xrange=[15,53],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='Modal Ages',margin=0.105,font_size=11, $ ; dimensions=[650,400]) ;ii = where(finite(modal_age_ind)) ;;p = plot(lat_ind[ii],tpause_theta_ind[ii],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) ;poly = polygon([mer_lat[215:-1],reverse(mer_lat[215:-1])],[tpm_theta_avg[215:-1]-tpm_theta_sd_avg[215:-1],reverse(tpm_theta_avg[215:-1])+reverse(tpm_theta_sd_avg[215:-1])],/fill_background,/data, $ ; linestyle=6,fill_color='orange',transparency=70) ;p = plot(mer_lat[215:-1],tpm_theta_avg[215:-1],color='dark orange',thick=3,/overplot) ;mc = alog10(modal_age_ind[ii])*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(lat_ind[ii],theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;tickvals = findgen(11)*0.25 ;tickname = strmid(strcompress(string(10^tickvals),/r),0,4) ;cb = COLORBAR(POSITION=[0.9015,0.23,0.9165,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10,title='Days') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_modal_age_ind.png' ;p = plot(indgen(2),/nodata,xrange=[13,55],yrange=[-30,150],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Tropopause Relative Theta (K)',title='SEAC$^4$RS Mean Age', $ ; margin=0.10,font_size=11,dimensions=[700,500]) ;p = plot([10,55],[0,0],linestyle=2,/overplot) ;ii = where(finite(mean_age_ind)) ;mc = mean_age_ind[ii]/365.*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(lat_ind[ii],theta_ind[ii]-tpause_theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;tickvals = findgen(13)*0.25 ;tickname = strmid(strcompress(string(tickvals),/r),0,4) ;cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='Years') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_mean_age_ind_tpause_rel.png' color_mult = 510. p = plot(indgen(2),/nodata,xrange=[15,53],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='NA Surface Source Fractions',margin=0.1, $ font_size=11,dimensions=[650,400]) ii = where(finite(mean_age_ind)) ;p = plot(lat_ind[ii],tpause_theta_ind[ii],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) poly = polygon([mer_lat[215:-1],reverse(mer_lat[215:-1])],[tpm_theta_avg[215:-1]-tpm_theta_sd_avg[215:-1],reverse(tpm_theta_avg[215:-1])+reverse(tpm_theta_sd_avg[215:-1])],/fill_background,/data, $ linestyle=6,fill_color='orange',transparency=70) p = plot(mer_lat[215:-1],tpm_theta_avg[215:-1],color='dark orange',thick=3,/overplot) mc = (1.-tr_source_frac[ii])*color_mult chk = where(mc gt 250,nchk) if nchk gt 0 then mc[chk] = 250 p = plot(lat_ind[ii],theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) tickvals = findgen(6)*0.1 tickname = strmid(strcompress(string(tickvals),/r),0,3) cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ title='Fraction') p.save,dir+'Plots/SEAC4RS/SEAC4RS_NH_frac_ind2.png' color_mult = 9. ;p = plot(indgen(2),/nodata,xrange=[15,53],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='Tropical Source Latitudes',margin=0.1,font_size=12, $ ; dimensions=[650,400]) ;for pp = 0, npeaks_source-1 do begin ; gd = where(tr_peak_ind eq pp and tr_source_frac gt 0.2 and co2_ind gt 0,ngd) ; if ngd gt 0 then begin ; for aa = 0, ngd-1 do begin ; mc = (tr_peaks[tr_peak_ind[gd[aa]]]+4)*color_mult ; chk = where(mc gt 250,nchk) ; if nchk gt 0 then mc[chk] = 250 ; chk = where(mc lt 0,nchk) ; if nchk gt 0 then mc[chk] = 0 ; p = plot(replicate(lat_ind[gd[aa]],2),replicate(theta_ind[gd[aa]],2),rgb_table=39,vert_colors=mc,symbol='o',sym_size=tr_source_frac[gd[aa]]*2,linestyle=6,/sym_filled,/overplot) ; endfor ; endif ;endfor ;tickvals = fix(tr_peaks)+4 ;tickname = strmid(strcompress(string(tickvals-4),/r),0,2) ;cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='Latitude') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_Tropical_source_lats.png' ; ;p = plot(indgen(2),/nodata,xrange=[15,53],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='NH Extratropical Source Latitudes',margin=0.1, $ ; font_size=12,dimensions=[650,400]) ;for pp = 0, npeaks_source-1 do begin ; gd = where(tr_peak_ind eq pp and tr_source_frac lt 0.8 and co2_ind gt 0,ngd) ; if ngd gt 0 then begin ; for aa = 0, ngd-1 do begin ; mc = (28 - (nh_peaks[nh_peak_ind[gd[aa]]]-nh_peaks[0]))*color_mult ; chk = where(mc gt 250,nchk) ; if nchk gt 0 then mc[chk] = 250 ; chk = where(mc lt 0,nchk) ; if nchk gt 0 then mc[chk] = 0 ; p = plot(replicate(lat_ind[gd[aa]],2),replicate(theta_ind[gd[aa]],2),rgb_table=39,vert_colors=mc,symbol='o',sym_size=(1.-tr_source_frac[gd[aa]])*2,linestyle=6,/sym_filled,/overplot) ; endfor ; endif ;endfor ;tickvals = fix(nh_peaks)-nh_peaks[0] ;tickname = strmid(strcompress(string(tickvals+nh_peaks[0]),/r),0,2) ;cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=reverse(rgb_39),BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='Latitude') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_NH_source_lats.png' color_mult = 0.1 ;i = 9 ;p = plot(indgen(2),/nodata,xrange=[13,55],yrange=[330,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)', $ ; title=strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4)+' Path Integrated Lifetime',margin=0.1,font_size=11,dimensions=[700,500]) ;ii = where(finite(reform(tau_ind[i,*]))) ;p = plot(lat_ind[ii],tpause_theta_ind[ii],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) ;mc = tau_ind[i,ii]*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(lat_ind[ii],theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;tickvals = findgen(11)*100. ;tickname = strmid(strcompress(string(tickvals),/r),0,4) ;cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='Years') color_mult = 125 ;p = plot(indgen(2),/nodata,xrange=[10,55],yrange=[310,505],xticklen=0.03,yticklen=0.03,xminor=1,xtitle='Latitude',ytitle='Theta (K)',title='SEAC$^4$RS Mean Age',margin=0.105,font_size=11, $ ; dimensions=[700,500]) ;ii = where(finite(mean_age_ind)) ;p = plot(lat_ind[ii],tpause_theta_ind[ii],color='orange',symbol='o',sym_size=0.5,linestyle=6,/overplot) ;mc = (alog10(mean_age_ind[ii])-1.5)*color_mult ;chk = where(mc gt 250,nchk) ;if nchk gt 0 then mc[chk] = 250 ;p = plot(lat_ind[ii],theta_ind[ii],rgb_table=39,vert_colors=mc,symbol='o',sym_size=1,linestyle=6,/sym_filled,/overplot) ;tickvals = (findgen(11)+1.5)/5.5 ;tickname = strmid(strcompress(string(tickvals),/r),0,4) ;cb = COLORBAR(POSITION=[0.91,0.23,0.925,0.77],tickname=tickname,tickvalues=tickvals*color_mult,RGB_TABLE=rgb_39,BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir,font_size=10, $ ; title='log10(Days)') n_levels = 30 ll = findgen(n_levels) lls = findgen(n_levels/4+1)*4 tickname = strmid(strcompress(string(fix(lls)),/r),0,2) z = 6 ;p = plot(indgen(2),/nodata,xrange=[0.5,2e3],/xlog,xtitle='Days',ytitle='Frequency',title='SEAC$^4$RS Transit Time Spectrum',dimensions=[700,500],margin=0.105,font_size=11) ;for z = 6, 18 do begin ;; p = plot(life_grid,gr_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;; p = plot(life_grid_linear,gr_linear_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;; p = plot(life_grid_linear,gr_linear_theta[*,0,z],color=colors[z-6],thick=1,/overplot) ;; for jj = 0, ni-1 do p = plot(life_grid_linear,gr_linear_theta[*,jj,z],color='grey',/overplot) ;; for jj = 0, ni-1 do p = plot(life_grid,gr_theta[*,jj,z],color='grey',/overplot) ;; p = plot(life_grid,gr_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;; for ki = 40, ng-1 do p = plot(life_grid_linear,g_linear_theta_all[*,ki,z],color='grey',/overplot) ;; p = plot(life_grid_linear,g_linear_theta_all[*,diff_g_min_i[z],z],color=colors[z-6],thick=2,/overplot) ;; p = plot(life_grid_linear,g_theta_best[*,z],color=colors[z-6],thick=2,/overplot) ; p = plot(life_grid_linear,g_best_source[*,z],color=colors[z-6],thick=2,/overplot) ;; p = plot(life_grid_linear,g_theta_best[*,z] + seas[*,0]*g_theta_best[*,z],color=colors[z-6],thick=2,/overplot) ; th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ; t = text(0.75,0.8-0.03*(z-6),th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors[z-6],font_size=11,/norm) ;endfor ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_theta1.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_theta.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_predef_theta.png' p = plot(indgen(2),/nodata,xrange=[0.5,1.5e2],yrange=[0,2.5],/xlog,xtitle='Days',ytitle='1e2*Frequency',title='NA Age Spectra',dimensions=[600,400],margin=[0.11,0.1,0.03,0.08],font_size=11) for z = 6, 18 do begin ; p = plot(life_grid_linear,g_best_source[*,z],color=colors[z-6],thick=3,/overplot) ; p = plot(life_grid_linear,g_best_source_tr[*,z],color=colors[z-6],thick=3,linestyle=2,/overplot) ; p = plot(life_grid_linear,1e2*g_best_source_nh[*,z],color=colors2[1.5*(z-5)],thick=2,/overplot) p = plot(life_grid_linear,1e2*g_best_source_nh[*,z],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ; p = plot(life_grid_linear,g_best_source_tr[*,z],color=colors2[1.5*(z-5)],thick=2,linestyle=2,/overplot) ; th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ; t = text(0.75,0.8-0.03*(z-6),th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors2[1.5*(z-5)],font_size=11,/norm) endfor p = plot(life_grid_linear,1e2*g_best_source_nh[*,10],color='deep sky blue',thick=2,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_nh_age_spectra_theta.png' ; ;p = plot(indgen(2),/nodata,xrange=[0,1.5],yrange=[0,2.8],xtitle='Years',ytitle='1e3*Frequency',title='Tropical Age Spectra',dimensions=[600,400],margin=[0.11,0.1,0.03,0.08],font_size=11) ;for z = 6, 18 do begin ;; p = plot(life_grid_linear,g_best_source[*,z],color=colors[z-6],thick=3,/overplot) ;; p = plot(life_grid_linear/365.,1e2*g_best_source_tr[*,z],color=colors2[1.5*(z-5)],thick=2,/overplot) ; p = plot(life_grid_linear/365.,1e3*g_best_source_tr[*,z],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;; p = plot(life_grid_linear,g_best_source_nh[*,z],color=colors[z-6],thick=3,linestyle=1,/overplot) ;; th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;; t = text(0.75,0.8-0.03*(z-6),th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors2[1.5*(z-5)],font_size=11,/norm) ;endfor ;p = plot(life_grid_linear/365.,1e3*g_best_source_tr[*,10],color='deep sky blue',thick=2,/overplot) ;;p = plot(life_grid_linear/365.,age_spec_tr_clams_seac4rs_interp[*,6],color='purple',thick=2,/overplot) ;tickvals = findgen(7)*20+theta_grid[6]-5 ;tickname = strmid(strcompress(string(tickvals),/r),0,3) ;cb = COLORBAR(POSITION=[0.78,0.4,0.8,0.82],tickname=tickname,tickvalues=1.5*(tickvals-theta_grid[6]+5),RGB_TABLE=rgb_3[70:-1,*],BORDER=1,ticklen=1,minor=0,ORIENTATION=1,/textpos,/tickdir, $ ; font_size=10,title='Theta (K)') ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_tr_age_spectra_theta.png' z = 9 th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;p = plot(indgen(2),/nodata,xrange=[0.5,2e2],xtitle='Days',ytitle='Frequency',title='SEAC$^4$RS NH Transit Time Spectrum',dimensions=[700,500],margin=0.11,font_size=11) ;p = plot(life_grid_linear,g_best_source_nh[*,z],color='blue',thick=2,/overplot) ;p = plot(life_grid_linear,age_spec_nh_clams_seac4rs_interp[*,z],color='red',thick=2,/overplot) ;p = plot(indgen(2),/nodata,xrange=[0,2],xtitle='Years',ytitle='Frequency',title='SEAC$^4$RS Tropical Transit Time Spectrum',dimensions=[700,500],margin=0.11,font_size=11) ;p = plot(life_grid_linear/365.,g_best_source_tr[*,z],color='blue',thick=2,/overplot) ;p = plot(life_grid_linear/365.,age_spec_tr_clams_seac4rs_interp[*,z],color='red',thick=2,/overplot) ;p = plot(indgen(2),/nodata,xrange=[0,10],xtitle='Years',/ylog,yrange=[1e-6,5e-3],ytitle='Frequency',title='SEAC$^4$RS Tropical Transit Time Spectrum',dimensions=[700,500],margin=0.11,font_size=11) ;p = plot(life_grid_linear/365.,g_best_source_tr[*,z],color='blue',thick=2,/overplot) ;p = plot(life_grid_linear/365.,age_spec_tr_clams_seac4rs_interp[*,z],color='red',thick=2,/overplot) ;p = plot(indgen(2),/nodata,xrange=[1,4e2],/xlog,xtitle='Days',yrange=[0,0.14],ytitle='Frequency*10',title='SEAC$^4$RS Age Spectra $\theta$='+th1+'K',dimensions=[700,500], $ ; margin=[0.095,0.08,0.03,0.07],font_size=11) ;for ki = ki_prof[z]-25, ki_prof[z]+20, 3 do p = plot(life_grid_linear,10.*g_linear_theta_all[*,ki,z],color='sky blue',/overplot) ;p = plot(life_grid_linear,10.*g_best_source_nh[*,z],color='medium orchid',thick=3,/overplot) ;p = plot(life_grid_linear,10.*g_best_source_tr[*,z],color='lime green',thick=3,/overplot) ;p = plot(life_grid_linear,10.*g_best_source[*,z],color='blue',thick=3,/overplot) ;;;for jj = 0, ni-2 do p = plot(life_grid,gr_theta[*,jj,z],color='medium orchid',thick=3,linestyle=2,/overplot) ;;for jj = 0, ni-2 do p = plot(life_grid,gr_theta[*,jj,z],color='dodger blue',thick=3,linestyle=2,/overplot) ;;;for jj = 0, 0 do p = plot(life_grid,gr_theta[*,jj,z],color='dodger blue',thick=3,/overplot) ;;p = plot(life_grid,gr_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;;t = text(0.2,0.76-0.03*(z-6),th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors[z-6],font_size=11,/norm) ;t = text(0.6,0.8,'Best total spectrum',color='blue',font_size=11,/norm) ;t = text(0.6,0.76,'From NH surface',color='medium orchid',font_size=11,/norm) ;t = text(0.6,0.72,'From Tropical surface',color='lime green',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectra_source_'+th1+'K.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_iter.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_iter_init.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_iter_init2.png' ;p = plot(indgen(2),/nodata,xrange=[0,35],yrange=[0,0.09],xticklen=0,yticklen=0,ytickfont_size=0,xtickfont_size=0,axis_style=1,dimensions=[400,400],margin=[0.095,0.08,0.03,0.07]) ;p = plot(life_grid_linear,10.*g_best_source[*,z],color='blue',thick=6,linestyle=2,/overplot) ;;p = plot(life_grid_linear,10.*g_best_source_tr[*,z],color='lime green',thick=6,/overplot) ;p = plot(life_grid_linear,10.*g_best_source_nh[*,z],color='medium orchid',thick=6,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_NA.png' ; ;p = plot(indgen(2),/nodata,xrange=[0,250],yrange=[0,0.09],xticklen=0,yticklen=0,ytickfont_size=0,xtickfont_size=0,axis_style=1,dimensions=[400,400],margin=[0.095,0.08,0.03,0.07]) ;p = plot(life_grid_linear,10.*g_best_source[*,z],color='blue',thick=6,linestyle=2,/overplot) ;;p = plot(life_grid_linear,10.*g_best_source_nh[*,z],color='medium orchid',thick=6,/overplot) ;p = plot(life_grid_linear,10.*g_best_source_tr[*,z],color='lime green',thick=6,/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_tropics.png' ;p = plot(indgen(2),/nodata,xrange=[0.1,1e4],/xlog,yrange=[0,0.04],xtitle='Days',ytitle='Frequency',title='SEAC$^4$RS Transit Time Spectrum',dimensions=[700,500],margin=0.095,font_size=11) ;p = plot(life_grid,gr_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;p = plot(365.*age_1,freq_1,color='red',thick=3,linestyle=2,/overplot) ;p = plot(365.*age_2,freq_2,color='orange',thick=3,linestyle=2,/overplot) ;t = text(0.15,0.81,th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors[z-6],font_size=11,/norm) ;t = text(0.15,0.76,'Hauck,JJA,370K,56N',color='red',font_size=11,/norm) ;t = text(0.15,0.72,'Hauck,SON,370K,56N',color='orange',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_Hauck_370K.png' ;z = 11 ;th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;p = plot(indgen(2),/nodata,xrange=[1,1e4],/xlog,yrange=[0,0.06],xtitle='Days',ytitle='Frequency',title='SEAC$^4$RS Transit Time Spectrum',dimensions=[700,500],margin=0.095,font_size=11) ;p = plot(life_grid,gr_theta[*,-1,z],color='blue',thick=3,/overplot) ;p = plot(365.*age_3,0.5*freq_3,color='red',thick=3,linestyle=2,/overplot) ;p = plot(365.*age_4,0.5*freq_4,color='orange',thick=3,linestyle=2,/overplot) ;p = plot(365.*age_5,2.5*freq_5,color='green',thick=3,linestyle=2,/overplot) ;p = plot(365.*age_6,2.5*freq_6,color='lime green',thick=3,linestyle=2,/overplot) ;t = text(0.15,0.81,th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color='blue',font_size=11,/norm) ;t = text(0.15,0.76,'Hauck,JJA,70hPa,55N',color='red',font_size=11,/norm) ;t = text(0.15,0.72,'Hauck,SON,70hPa,55N',color='orange',font_size=11,/norm) ;t = text(0.15,0.67,'Ploeger,JJA,400K,60N',color='green',font_size=11,/norm) ;t = text(0.15,0.63,'Ploeger,SON,400K,60N',color='lime green',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_Hauck_Ploeger_400K_55N.png' ; ;p = plot(indgen(2),/nodata,xrange=[1,1e4],/xlog,yrange=[0,0.06],xtitle='Days',ytitle='Frequency',title='SEAC$^4$RS Transit Time Spectrum',dimensions=[700,500],margin=0.095,font_size=11) ;p = plot(life_grid,gr_theta[*,-1,z],color='blue',thick=3,/overplot) ;p = plot(365.*age_5,2.5*freq_5,color='red',thick=3,linestyle=2,/overplot) ;p = plot(365.*age_6,2.5*freq_6,color='orange',thick=3,linestyle=2,/overplot) ;t = text(0.15,0.81,th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color='blue',font_size=11,/norm) ;t = text(0.15,0.76,'Ploeger,JJA,400K,60N',color='red',font_size=11,/norm) ;t = text(0.15,0.72,'Ploeger,SON,400K,60N',color='orange',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectrum_'+th1+'K_Ploeger_400K.png' ;p = plot(indgen(2),/nodata,xrange=[0.1,5e3],/xlog,yrange=[0,0.11],xtitle='Days',ytitle='Frequency',title='SEAC4RS Transit Spectrum') ;;for jj = 0, nii-1 do p = plot(life_grid,gr_ind[*,jj,aa],color='grey',/overplot) ;p = plot(life_grid,gr_ind[*,-1,aa],color='purple',thick=2,/overplot) z = 11 ;p = plot(indgen(2),/nodata,xrange=[1,3e2],/xlog,xtitle='Days',ytitle='1e2*Frequency',title='Age Spectra 400-410 K',dimensions=[700,500],margin=0.095,font_size=11) ;zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;for n = 0, nn-1 do p = plot(life_grid_linear,1e2*g_best_ind[*,zi_th[n]],color='sky blue',thick=2,/overplot) ;for n = 0, nn-1 do p = plot(life_grid_linear,1e2*g_best_ind_nh[*,zi_th[n]],color='medium orchid',thick=2,/overplot) ;for n = 0, nn-1 do p = plot(life_grid_linear,1e2*g_best_ind_tr[*,zi_th[n]],color='lime green',thick=2,/overplot) ;p = plot(life_grid_linear,1e2*g_best_source[*,z],color='blue',thick=2,/overplot) ;t = text(0.6,0.82,'Theta average spectrum',color='blue',font_size=11,/norm) ;t = text(0.6,0.78,'Individual total spectra',color='sky blue',font_size=11,/norm) ;t = text(0.6,0.74,'From NH surface',color='medium orchid',font_size=11,/norm) ;t = text(0.6,0.7,'From Tropical surface',color='lime green',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_age_spectra_ind.png' ;p = plot(indgen(2),/nodata,xrange=[0,10],yrange=[0,0.2],ytitle='Frequency',title='SEAC4RS Transit Spectrum') ;;for i = 0, ng-1 do p = plot(life_grid/365.,g[*,0,i],color='green',thick=2,/overplot) ;for i = 0, ng-1 do p = plot(life_grid/365.,g_theta_all[*,i,6],color='orange',thick=2,/overplot) ;p = plot(indgen(2),/nodata,xrange=[1,1e5],/xlog,yrange=[0,1.05],xtitle='Lifetime ($\tau$, Days)',ytitle='Boundary Layer Fraction ($\mu$)',title='$\mu$ vs $\tau$',dimensions=[600,400], $ ; margin=[0.11,0.1,0.03,0.08],font_size=11) ;;p = plot([1e3,1e5],[1,1],linestyle=2,/overplot) ;for z = 6, 18 do begin ;; for jj = 0, 1 do p = plot(ci_lifetimes_linear_theta_profile[jj,1:-1,z],trop_fracs_linear_theta_gradj[jj,z,1:-1],color='lime green',symbol='o',/sym_filled,linestyle=6,/overplot) ;; for jj = 0, ni-1 do p = plot(ci_lifetimes_linear_theta_profile[jj,1:-1,z],trop_fracs_linear_theta_gradj[jj,z,1:-1],color='grey',symbol='s',linestyle=6,/overplot) ;; p = plot(ci_lifetimes_linear_theta_profile[-1,1:-1,z],trop_fracs_linear_theta_gradj[-1,z,1:-1],color=colors[z-6],symbol='s',/sym_filled,linestyle=6,/overplot) ;; p = plot(tau_profile_source[1:-1,z],tf_best_source[z,1:-1],color=colors2[1.5*(z-5)],symbol='o',/sym_filled,linestyle=6,/overplot) ;; for ki = 0, 50 do p = plot(ci_lifetimes_theta_profile[-1,1:-1,z],bl_frac_gr_theta_region[ki,0,z,1:-1],color='grey',symbol='o',linestyle=6,/overplot) ;; for jj = 0, 1 do p = plot(life_grid_linear,trop_fracs_linear_theta_fit[*,jj,z],color='grey',/overplot) ;; p = plot(life_grid,trop_frac_gr_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ;; p = plot(life_grid_linear,trop_frac_gr_linear_theta[*,-1,z],color=colors[z-6],thick=3,/overplot) ; p = plot(life_grid_linear,tfk_best_source[*,z],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;; p = plot(life_grid_linear,trop_frac_gr_linear_theta[*,-1,z-1],color=colors[z-1-6],thick=3,linestyle=2,/overplot) ;; for jj = 0, 1 do p = plot(life_grid_linear,trop_frac_gr_linear_theta[*,jj,z],color='lime green',linestyle=2,/overplot) ;; for jj = 0, ni-1 do p = plot(life_grid_linear,trop_frac_gr_linear_theta[*,jj,z],color='grey',thick=2,/overplot) ;; for jj = 0, ni-1 do p = plot(life_grid,trop_frac_gr_theta[*,jj,z],color='grey',linestyle=2,/overplot) ;; for ki = 0, ng-1 do p = plot(life_grid_linear,trop_frac_g_linear_theta[*,ki,z],color='grey',/overplot) ;; p = plot(life_grid_linear,trop_frac_g_linear_theta[*,diff_g_min_i[z],z],color='medium orchid',thick=2,/overplot) ;; th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;; t = text(0.75,0.4-0.03*(z-12),th1+'-'+strcompress(string(fix(theta_grid[z]+dth/2)),/r)+'K',color=colors2[1.5*(z-5)],font_size=11,/norm) ;endfor ;p = plot(life_grid_linear,tfk_best_source[*,10],color='deep sky blue',thick=2,/overplot) ;for z = 6, 18 do for i = 1, nciv-1 do s = symbol(tau_profile_source[i,z],tf_best_source[z,i],'o',sym_color=colors2[i],/sym_filled,/data) ;;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_Trop_frac_vs_lifetime_linear.png' ;;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_Trop_frac_vs_lifetime_linear2.png' ;;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_linear1.png' ;;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_linear.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source.png' ;p = plot(indgen(2),/nodata,xrange=[2e3,1e5],/xlog,yrange=[0.9,1.002],xtitle='Lifetime (Days)',ytitle='Boundary Layer Fraction ($\mu$)',dimensions=[550,500], $ ; margin=[0.13,0.09,0.02,0.02],font_size=12) ;for z = 6, 18 do p = plot(life_grid_linear,tfk_best_source[*,z],rgb_table=3,vert_colors=1.4*(theta_grid[z]-theta_grid[6])+70,thick=2,/overplot) ;for z = 6, 18 do for i = 1, nciv-1 do s = symbol(tau_profile_source[i,z],tf_best_source[z,i],'o',sym_color=colors2[i],/sym_filled,/data) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_inset.png' z = 8 th1 = strcompress(string(fix(theta_grid[z]-dth/2)),/r) ;p = plot(indgen(2),/nodata,xrange=[0.5,1e5],/xlog,yrange=[0,1.05],xtitle='Lifetime (Days)',ytitle='Boundary Layer Fraction ($\mu$)',title='SEAC$^4$RS $\mu$ vs. $\tau$ $\theta$='+th1+'K', $ ; dimensions=[700,500],margin=0.08,font_size=11) ;p = plot([1e3,1e5],[1,1],linestyle=2,/overplot) ;;;p = plot(ci_lifetimes_linear_theta_profile[0,1:-1,z],trop_fracs_linear_theta_gradj[0,z,1:-1],color=colors[z-6],symbol='s',/sym_filled,linestyle=6,/overplot) ;;;p = plot(life_grid_linear,trop_frac_gr_linear_theta[*,0,z],color=colors[z-6],thick=3,/overplot) ;;;for ki = 0, ng-1 do p = plot(ci_lifetimes_linear_theta_profile[0,1:-1,z],tf_gradj[ki,z,1:-1],color='grey',symbol='s',linestyle=6,/overplot) ;;;for ki = 100, 120 do p = plot(tau_k[ki,z,1:-1],tf_gradj[ki,z,1:-1],color='grey',symbol='s',linestyle=6,/overplot) ;;;p = plot(tau_k[20,z,1:-1],tf_gradj[1,z,1:-1],color='green',symbol='s',/sym_filled,linestyle=6,/overplot) ;;;p = plot(tau_k[-1,z,1:-1],tf_gradj[-1,z,1:-1],color='red',symbol='s',/sym_filled,linestyle=6,/overplot) ;;;p = plot(tau_profile[1:-1,z],tf_gradj_best[z,1:-1],color='blue',symbol='s',/sym_filled,linestyle=6,/overplot) ;;for ki = ki_prof[z]-15, ki_prof[z]+10, 2 do p = plot(life_grid_linear,tfk[*,ki,z],color='sky blue',/overplot) ;for ki = ki_prof[z]-15, ki_prof[z]+10, 2 do p = plot(life_grid_linear,tfk[*,ki,z],color='sky blue',thick=2,/overplot) ;;p = plot(life_grid_linear,tf_clams[*,z],thick=1,color='red',/overplot) ;;for pp = 0, npeaks_source-1 do for pp2 = 0, npeaks_source-1 do $ ;; p = plot(tau_profile_source[1:-1,z],tf_clams_all[pp,0,pp2,0,z,1:-1],color='red',symbol='o',/sym_filled,sym_size=0.5,linestyle=6,/overplot) ;;p = plot(replicate(tau_profile_source[8,z],npeaks_source^2),tf_clams_all[*,0,*,0,z,8],color='purple',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;;p = plot(tau_profile_source[1:-1,z],tf_clams_all[tr_peak[z],tr_width[z],nh_peak[z],nh_width[z],z,1:-1],color='orange',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,z],thick=2,color='blue',/overplot) ;for ki = ki_prof[z], ki_prof[z] do for pp = 0, npeaks_source-1 do for pp2 = 0, npeaks_source-1 do $ ; p = plot(tau_profile_source[1:-1,z-1],tf_surf[pp,0,pp2,0,ki,z,1:-1],symbol='o',/sym_filled,color='orange',sym_size=0.5,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z-1],tf_best_source[z,1:-1],color='blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z],tf_best_source[z,1:-1],color='red',symbol='s',/sym_filled,linestyle=6,/overplot) ;for i = 1, nciv-1 do p = plot(tau_profile_source[i,z-1:z],replicate(tf_best_source[z,i],2),color='red',/overplot) ;;;for ki = 0, ng-1 do p = plot(life_grid_linear,tfk_fit[*,ki,z],color='grey',/overplot) ;;;p = plot(life_grid_linear,tfk_fit[*,20,z],color='green',/overplot) ;;;p = plot(life_grid_linear,tfk_fit[*,-1,z],color='red',/overplot) ;;;p = plot(life_grid_linear,tfk_best[*,z],color='medium orchid',thick=2,/overplot) ;;t = text(0.17,0.8,'$\mu_G$',color='sky blue',font_size=12,/norm) ;;t = text(0.17,0.76,'$\mu_i$',color='orange',font_size=12,/norm) ;;t = text(0.17,0.72,'Best fit',color='blue',font_size=12,/norm) ;;t = text(0.17,0.68,'Adjusted $\tau_i$',color='red',font_size=12,/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_'+th1+'K.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_'+th1+'K_1.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_'+th1+'K_2.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_'+th1+'K_3.png' ; ;p = plot(indgen(2),/nodata,xrange=[2e3,7e4],/xlog,yrange=[0.91,1.005],xtitle='Lifetime (Days)',ytitle='Boundary Layer Fraction ($\mu$)',dimensions=[550,500], $ ; margin=[0.13,0.09,0.02,0.02],font_size=12) ;p = plot([1e3,1e5],[1,1],linestyle=2,/overplot) ;for ki = ki_prof[z]-25, ki_prof[z]+10, 3 do p = plot(life_grid_linear,tfk[*,ki,z],color='sky blue',thick=2,/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,z],thick=3,color='blue',/overplot) ;for ki = ki_prof[z], ki_prof[z] do for pp = 0, npeaks_source-1 do for pp2 = 0, npeaks_source-1 do $ ; p = plot(tau_profile_source[1:-1,z-1],tf_surf[pp,0,pp2,0,ki,z,1:-1],symbol='o',/sym_filled,color='orange',sym_size=1,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z-1],tf_best_source[z,1:-1],color='blue',symbol='o',/sym_filled,sym_size=1,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z],tf_best_source[z,1:-1],color='red',symbol='s',/sym_filled,sym_size=2,linestyle=6,/overplot) ;for i = 1, nciv-1 do p = plot(tau_profile_source[i,z-1:z],replicate(tf_best_source[z,i],2),color='red',/overplot) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_inset_'+th1+'K.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_source_inset_'+th1+'K_1.png' ;p = plot(indgen(2),/nodata,xrange=[1,1e5],/xlog,yrange=[1e-2,7],/ylog,xtitle='Lifetime (Days)',ytitle='Boundary Layer Fraction ($\mu$)',title='SEAC$^4$RS $\mu$ vs. $\tau$ $\theta$='+th1+'K', $ ; dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([1,1e5],[1,1],linestyle=2,/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,z],thick=2,color='blue',/overplot) ;p = plot(tau_profile_source[1:-1,z],tf_best_source[z,1:-1],color='red',symbol='s',/sym_filled,linestyle=6,/overplot) ;zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;;if nn gt 0 then for i = 1, 20 do p = plot(replicate(tau_profile_source[i,z],nn),tf_ind_init[zi_th,i],color=colors2[i],symbol='o',sym_size=0.5,linestyle=6,/overplot) ;if nn gt 0 then for i = 1, 20 do p = plot(replicate(tau_profile_source[i,z],nn),tf_ind_adj[zi_th,i],color=colors2[i],symbol='o',sym_size=0.5,linestyle=6,/overplot) ;gd = where(finite(tf_ind_init[zi_th,3])) ;si = sort(tf_ind_init[zi_th[gd],3]) ;p = plot(tau_profile_source[1:-1,z],reform(tf_ind_init[zi_th[gd[si[0]]],1:-1]),color='medium orchid',symbol='o',/sym_filled,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z],reform(tf_ind_adj[zi_th[gd[si[0]]],1:-1]),color=colors2[i],symbol='tu',sym_size=1,/sym_filled,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z],tf_ind_adj_mean[z,1:-1],color='purple',symbol='s',/sym_filled,linestyle=6,/overplot) gd = where(finite(tf_ind_init[*,3])) si = reverse(sort(tf_ind_init[gd,3])) aa = gd[si[6]] zii = round(zi[aa]) ;p = plot(indgen(2),/nodata,xrange=[1,1.2e3],/xlog,yrange=[2e-2,4],/ylog,xtitle='Lifetime ($\tau$, Days)',ytitle='Boundary Layer Fraction ($\mu$)', $ ; title='Single Location $\mu$ vs. $\tau$',dimensions=[700,500],margin=0.09,font_size=11) ;poly = polygon([0.4,2e3,2e3,0.4],[1,1,5,5],/fill_background,/data,linestyle=6,fill_color='orange',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot([0.4,1e5],[1,1],linestyle=2,/overplot) ;low = where(0.2*life_grid_linear lt 1) ;lg = 0.25*life_grid_linear ;lg[low] = 1.0 ;poly = polygon([0.4,lg,0.4],[2e-2,tfk_best_source[*,zii],1],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;lg = 4.0*life_grid_linear ;low = where(lg lt 2e3) ;poly = polygon([lg[low],2e3,2e3],[tfk_best_source[low,zii],1,2e-2],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot([0.4,1e5],[1,1],linestyle=2,/overplot) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_ind_init[aa,i],sym_color=colors4[i],'o',sym_size=2,sym_thick=2,/data) ;for i = 1, 10 do t = text(0.71,0.42-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],/norm) ;p = plot(life_grid_linear,tfk_best_ind[*,aa],thick=2,color='blue',/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,zii],thick=2,color='blue',linestyle=1,/overplot) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_best_ind[aa,i],sym_color=colors4[i],'s',sym_size=1.5,/data) ;for i = 1, 10 do p = symbol(tau_ind[i,aa],tf_best_ind[aa,i],sym_color=colors4[i],'s',sym_size=1.5,/sym_filled,/data) ;p = symbol(0.51,0.21,'s',sym_size=1.5,/norm) ;t = text(0.53,0.2,'optimized $\mu$',font_size=11,/norm) ;p = symbol(0.51,0.17,'s',sym_size=1.5,/sym_filled,/norm) ;t = text(0.53,0.16,'optimized $\tau$',font_size=11,/norm) ;t = text(0.12,0.81,'Trace gases with',font_size=11,color='orange',/norm) ;t = text(0.12,0.77,'polluted BL source',font_size=11,color='orange',/norm) ;t = text(0.12,0.645,'$\tau$ beyond range',font_size=11,color='lime green',/norm) ;t = text(0.12,0.605,'for theta layer',font_size=11,color='lime green',/norm) ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers_short1.png' ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers_short2.png' ;;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers_short3.png' ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers_short4.png' ;p = plot(indgen(2),/nodata,xrange=[1,1.2e3],/xlog,yrange=[2e-2,4],/ylog,xtitle='Lifetime ($\tau$, Days)',ytitle='Boundary Layer Fraction ($\mu$)', $ ; title='Single Location $\mu$ vs. $\tau$',dimensions=[700,500],margin=0.09,font_size=11) ;poly = polygon([0.4,2e3,2e3,0.4],[1,1,5,5],/fill_background,/data,linestyle=6,fill_color='orange',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot([0.4,1e5],[1,1],linestyle=2,/overplot) ;low = where(0.2*life_grid_linear lt 1) ;lg = 0.25*life_grid_linear ;lg[low] = 1.0 ;poly = polygon([0.4,lg,0.4],[2e-2,tfk_best_source[*,zii],1],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;lg = 4.0*life_grid_linear ;low = where(lg lt 2e3) ;poly = polygon([lg[low],2e3,2e3],[tfk_best_source[low,zii],1,2e-2],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;;p = plot(life_grid_linear,tfk_best_ind_init[*,aa],thick=2,color='blue',linestyle=2,/overplot) ;kii = where(diffs_ind_tot[tr_peak_ind[aa],nh_peak_ind[aa],*,aa] lt 1) ;;p = plot(life_grid_linear,tfk[*,kii[0],round(zii)],thick=2,color='sky blue',/overplot) ;;p = plot(life_grid_linear,tfk[*,kii[-1],round(zii)],thick=2,color='sky blue',/overplot) ;poly = polygon([life_grid_linear,reverse(life_grid_linear)],[tfk[*,kii[0],round(zii)],reverse(tfk[*,kii[-1],round(zii)])],/fill_background,/data,linestyle=6,fill_color='sky blue',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot(life_grid_linear,tfk_best_ind[*,aa],thick=2,color='blue',/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,zii],thick=2,color='blue',linestyle=1,/overplot) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_ind_init[aa,i],sym_color=colors4[i],'o',sym_size=2,sym_thick=2,/data) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_ind_adj[aa,i],sym_color=colors4[i],'tu',sym_size=1.5,sym_thick=2,/data) ;;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],scale_tf[i,aa]*tf_ind_adj[aa,i],sym_color=colors4[i],'tl',sym_size=1.5,sym_thick=2,/data) ;;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_ind_all_scale[tr_peak_ind[aa],nh_peak_ind[aa],aa,i],sym_color=colors4[i],'td',sym_size=1.5,sym_thick=2,/data) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],scale_tf_tau[i,aa]*scale_tf[i,aa]*tf_ind_adj[aa,i],sym_color=colors4[i],'td',sym_size=1.5,sym_thick=2,/data) ;;for i = 1, 10 do p = symbol(tau_ind_init[i,aa],tf_best_ind_init[aa,i],sym_color=colors4[i],'tu',sym_size=1.5,/sym_filled,/data) ;for i = 1, 10 do p = symbol(tau_profile_source[i,zii],tf_best_ind[aa,i],sym_color=colors4[i],'s',sym_size=1.5,/data) ;for i = 1, 10 do p = symbol(tau_ind[i,aa],tf_best_ind[aa,i],sym_color=colors4[i],'s',sym_size=1.5,/sym_filled,/data) ;for i = 1, 10 do t = text(0.71,0.42-(i-1)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i],/norm) ;p = symbol(0.51,0.37,'o',sym_size=2,sym_thick=2,/norm) ;t = text(0.53,0.36,'Initial',font_size=11,/norm) ;p = symbol(0.51,0.33,'tu',sym_size=1.5,sym_thick=2,/norm) ;t = text(0.53,0.32,'with $\itS_{inorm}$',font_size=11,/norm) ;p = symbol(0.51,0.29,'tl',sym_size=1.5,sym_thick=2,/norm) ;t = text(0.53,0.28,'and $\itS_{i\mu}$',font_size=11,/norm) ;p = symbol(0.51,0.25,'td',sym_size=1.5,sym_thick=2,/norm) ;t = text(0.53,0.24,'and $\itS_{i\tau}$',font_size=11,/norm) ;p = symbol(0.51,0.21,'s',sym_size=1.5,/norm) ;t = text(0.53,0.2,'optimized $\mu$',font_size=11,/norm) ;p = symbol(0.51,0.17,'s',sym_size=1.5,/sym_filled,/norm) ;t = text(0.53,0.16,'optimized $\tau$',font_size=11,/norm) ;t = text(0.12,0.81,'Trace gases with',font_size=11,color='orange',/norm) ;t = text(0.12,0.77,'polluted BL source',font_size=11,color='orange',/norm) ;t = text(0.12,0.64,'$\tau$ beyond range',font_size=11,color='lime green',/norm) ;t = text(0.12,0.6,'for theta layer',font_size=11,color='lime green',/norm) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers.png' ;print,reform(diffs_ind_tot[tr_peak_ind[aa]-1,nh_peak_ind[aa]-1,*,aa]),tr_peak_ind[aa],nh_peak_ind[aa],ki_ind[aa],mean_age_ind[aa],modal_age_ind[aa],lat_ind[aa],theta_ind[aa],co2_ind[aa] ;print,tau_ind[1:10,aa],reform(tf_best_ind[aa,1:10]),reform(scale_nh[aa,1:10]),scale_tf_tau[1:10,aa] ;p = plot(indgen(2),/nodata,xrange=[2e3,1e5],/xlog,yrange=[0.95,1.015],xtitle='Lifetime ($\tau$, Days)',ytitle='Boundary Layer Fraction ($\mu$)', $ ; title='Single Location $\mu$ vs. $\tau$',dimensions=[700,500],margin=0.095,font_size=11) ; poly = polygon([8e2,1e5,1e5,8e2],[1,1,5,5],/fill_background,/data,linestyle=6,fill_color='orange',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot([0.4,1e5],[1,1],linestyle=2,/overplot) ;low = where(0.2*life_grid_linear lt 1) ;lg = 0.25*life_grid_linear ;lg[low] = 1.0 ;poly = polygon([0.4,lg,7e4,0.4],[2e-2,tfk_best_source[*,zii],1,1],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;lg = 4.0*life_grid_linear ;low = where(lg lt 1e5) ;poly = polygon([lg[low],1e5,1e5],[tfk_best_source[low,zii],1,2e-2],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=2) ;kii = where(diffs_ind_tot[tr_peak_ind[aa],nh_peak_ind[aa],*,aa] lt 1) ;poly = polygon([life_grid_linear,reverse(life_grid_linear)],[tfk[*,kii[0],round(zii)],reverse(tfk[*,kii[-1],round(zii)])],/fill_background,/data,linestyle=6,fill_color='sky blue',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot(life_grid_linear,tfk_best_ind[*,aa],thick=2,color='blue',/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,zii],thick=2,color='blue',linestyle=1,/overplot) ;for i = 11, 20 do begin ; tf = tf_ind_init[aa,i] ; if i eq 13 then tf = 1.012 ; p = symbol(tau_profile_source[i,zii],tf,sym_color=colors4[i-10],'o',sym_size=2,sym_thick=2,/data) ;endfor ;;for i = 11, 20 do begin ;; tf = tf_ind_adj[aa,i] ;; if i eq 13 then tf = 1.015 ;; p = symbol(tau_profile_source[i,zii],tf,sym_color=colors4[i-10],'tu',sym_size=1.5,sym_thick=2,/data) ;;endfor ;for i = 11, 20 do p = symbol(tau_profile_source[i,zii],scale_tf[i,aa]*tf_ind_adj[aa,i],sym_color=colors4[i-10],'tl',sym_size=1.5,sym_thick=2,/data) ;for i = 11, 20 do p = symbol(tau_profile_source[i,zii],scale_tf_tau[i,aa]*scale_tf[i,aa]*tf_ind_adj[aa,i],sym_color=colors4[i-10],'td',sym_size=1.5,sym_thick=2,/data) ;for i = 11, 20 do p = symbol(tau_profile_source[i,zii],tf_best_ind[aa,i],sym_color=colors4[i-10],'s',sym_size=1.5,/data) ;for i = 11, 20 do p = symbol(tau_ind[i,aa],tf_best_ind[aa,i],sym_color=colors4[i-10],'s',sym_size=1.5,/sym_filled,/data) ;for i = 11, 20 do t = text(0.72,0.43-(i-10)*0.03,strmid(ci_vars_dc8[i],0,strlen(ci_vars_dc8[i])-4),color=colors4[i-10],/norm) ;t = text(tau_profile_source[13,zii],1.008,'(1.06)',color=colors4[3],/data,alignment=0.5) ;p.save,dir+'Plots/SEAC4RS/SEAC4RS_BL_frac_vs_lifetime_outliers2.png' ;p = plot(indgen(2),/nodata,xrange=[1,1e5],/xlog,yrange=[1.1,5],xtickfont_size=0,title='$\mu$ vs. $\tau$',dimensions=[700,100],margin=0.09,font_size=11) ;for i = 1, 20 do p = symbol(tau_profile_source[i,zii],tf_ind_init[aa,i],sym_color=colors2[i],'o',sym_size=1.5,/data) ;p = plot(indgen(2),/nodata,xrange=[1,1e5],/xlog,yrange=[0,10],xtitle='Lifetime ($\tau$, Days)',ytitle='Boundary Layer Fraction ($\mu$)',title='SEAC$^4$RS $\mu$ vs. $\tau$', $ ; dimensions=[700,500],margin=0.09,font_size=11) ;p = plot([1,1e5],[1,1],linestyle=2,/overplot) ;p = plot(life_grid_linear,tfk_best_source[*,zi[aa]],thick=2,color='blue',/overplot) ;p = plot(tau_profile_source[1:-1,zi[aa]],tf_best_source[zi[aa],1:-1],color='red',symbol='s',/sym_filled,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,zi[aa]],reform(tf_ind_init[aa,1:-1]),color='orange',symbol='o',linestyle=6,sym_size=1.5,/overplot) ;p = plot(tau_profile_source[1:-1,zi[aa]],reform(tf_ind_adj[aa,1:-1]),color='medium orchid',symbol='o',/sym_filled,linestyle=6,/overplot) ;p = plot(tau_ind[1:-1,aa],reform(tf_best_ind[aa,1:-1]),color='purple',symbol='s',/sym_filled,linestyle=6,/overplot) ;p = plot(life_grid_linear,tfk_best_ind[*,aa],thick=2,color='sky blue',/overplot) ;;zi_th = where(theta_ind ge theta_grid[z]-dth/2. and theta_ind lt theta_grid[z]+dth/2.,nn) ;;if nn gt 0 then for i = 1, 20 do p = plot(replicate(tau_profile_source[i,z],nn),tf_ind_init[zi_th,i],color=colors2[i],symbol='o',sym_size=0.5,linestyle=6,/overplot) ;gd = where(finite(tf_ind_init[zi_th,3])) ;si = sort(tf_ind_init[zi_th[gd],3]) ;p = plot(tau_profile_source[1:-1,z],reform(tf_ind_init[zi_th[gd[si[-1]]],1:-1]),color='medium orchid',symbol='o',/sym_filled,linestyle=6,/overplot) ;p = plot(tau_profile_source[1:-1,z],reform(tf_ind_adj[zi_th[gd[si[-1]]],1:-1]),color=colors2[i],symbol='tu',sym_size=1,/sym_filled,linestyle=6,/overplot) ;tau = fltarr(nciv) ;for i = 1, 20 do tau[i] = interpolate(tau_profile_source[i,*],zi[aa]) ;p = plot(indgen(2),/nodata,xrange=[1,1e5],/xlog,/ylog,yrange=[5e-2,50]) ;p = plot([1,1e5],[1,1],linestyle=2,/overplot) ;for i = 1, 20 do p = symbol(tau[i],tau_ind[i,aa]/tau[i],sym_color=colors2[i],'o',/sym_filled,sym_size=1,/data) ;p = plot(indgen(2),/nodata,xrange=[2012,2014],xtitle='Year',yrange=[380,408],ytitle='ppmv',title='GML Surface CO$_2$ Time Series',dimensions=[700,500],margin=0.09,font_size=11) ;;p = plot(lef_time_raw,lef_co2_raw,thick=2,color=colors[0],linestyle=1,/overplot) ;;p = plot(lef_time,co2_lef_mon,thick=2,color=colors[0],linestyle=2,/overplot) ;;p = plot(co2_time,co2_lef,thick=2,color=colors[5],/overplot) ;;for y = 0, 5 do p = plot(co2_time_ext,co2_surface_flask[*,y],thick=3,color=colors[y+1],/overplot) ;for y = 60, 84, 8 do p = plot(co2_time_ext,co2_surface_interp[*,y],thick=3,color=colors3[(y-52)/8],/overplot) ;for y = 1, 7, 2 do t = text(0.4,0.45-y/2*0.04,strcompress(string(fix(nh_peaks[y])),/r)+'$\deg$N',color=colors3[y],font_size=11,alignment=0.5,/norm) ;;p = plot(co2_time_ext,co2_ct_glb_zm[*,13],thick=3,color='medium orchid',/overplot) ;;for pp = 0, 4, 2 do p = plot(co2_time_ext,co2_surf[*,pp,0,1],thick=3,color='orange',/overplot) ;;for pp = 0, 4, 2 do p = plot(co2_time_ext,co2_surf[*,pp,0,0],thick=3,color='gold',/overplot) ;;p = plot(co2_time_ext,co2_ct_nam_zm[*,13],thick=3,color='pink',/overplot) ;;p = plot(co2_time_ext,co2_ct_combo_zm[*,13],thick=3,color='pink',/overplot) ;;p = plot(co2_time,co2_sgp,thick=2,color=colors[4],/overplot) ;;p = plot(co2_time,co2_key,thick=2,color=colors[3],/overplot) ;;p = plot(co2_time,co2_mlo,thick=2,color=colors[2],/overplot) ;;p = plot(co2_time,co2_brw,thick=2,color=colors[6],/overplot) ;;p = plot(co2_time,co2_smo,thick=2,color=colors[1],/overplot) ;poly = polygon([min(er2_yrfrac),max(er2_yrfrac),max(er2_yrfrac),min(er2_yrfrac)],[380,380,408,408],/fill_background,/data,linestyle=6,fill_color='purple',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ;p = plot(replicate(er2_yrfrac_mean,2),[380,410],thick=2,linestyle=2,color='purple',/overplot) ;;t = text(0.17,0.32,'BRW',color=colors[6],font_size=11,/norm) ;t = text(0.17,0.29,'LEF',color=colors[5],font_size=11,/norm) ;t = text(0.17,0.26,'SGP',color=colors[4],font_size=11,/norm) ;t = text(0.17,0.23,'KEY',color=colors[3],font_size=11,/norm) ;t = text(0.17,0.2,'MLO',color=colors[2],font_size=11,/norm) ;t = text(0.17,0.17,'SMO',color=colors[1],font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/CO2_surface_time_series.png' ;p = plot(indgen(2),/nodata,xrange=[2013,2014],xtickname=months,xminor=0,xmajor=13,xtitle='2013',yrange=[388,406],ytitle='ppm',dimensions=[750,400],margin=[0.09,0.1,0.02,0.09],font_size=12, $ ; title='CarbonTracker NAM BL CO$_2$ Time Series') ;for y = 60, 84, 4 do p = plot(co2_time_ext,co2_surface_interp[*,y],thick=3,color=colors3[(y-60)/4+1],/overplot) ;for y = 1, 7 do t = text(0.4,0.55-y*0.05,strcompress(string(fix(nh_peaks[y])),/r)+'$\deg$N',color=colors3[y],font_size=12,alignment=0.5,/norm) ;chk = where(co2_ind gt 0) ;poly = polygon([er2_yrfrac[chk[0]],max(er2_yrfrac),max(er2_yrfrac),er2_yrfrac[chk[0]]],[380,380,408,408],/fill_background,/data,linestyle=6,fill_color='purple',transparency=80, $ ; PATTERN_ORIENTATION=45,pattern_spacing=2) ;t = text(0.67,0.8,'SEAC$^4$RS',color='purple',font_size=11,/norm) ;t = text(0.67,0.76,'w/ER-2 CO$_2$',color='purple',font_size=11,/norm) ;p.save,dir+'Plots/SEAC4RS/CO2_CarbonTracker_NAM_BL_time_series.png' slats = [-14.3,19.5,25.8,46] slons = [187,205,280,270] alats = [-14.3,13,41,53] alons = [190,301,236,350] ;map = MAP('Equirectangular',label_position=0,linestyle='dotted',title='Measurement Locations',limit=[-25,170,60,360],center_longitude=260,/box_axes,margin=0.1,dimensions=[700,500]) ;m1 = MAPCONTINENTS(fill_color='white') ;m1 = MAPCONTINENTS(/usa) ;poly = polygon([248,275,275,248],[25,25,60,60],/fill_background,/data,linestyle=6,fill_color='medium orchid',transparency=70) ;poly = polygon([170,360,360,170],[25,25,60,60],/fill_background,/data,linestyle=6,fill_color='medium orchid',transparency=80,PATTERN_ORIENTATION=45,pattern_spacing=3) ;poly = polygon([170,360,360,170],[-25,-25,25,25],/fill_background,/data,linestyle=6,fill_color='lime green',transparency=70) ;for y = 0, 3 do s = symbol(slons[y],slats[y],/DATA,'star',sym_color='red',sym_size=2.5,/sym_filled) ;for y = 0, 3 do s = symbol(alons[y],alats[y],/DATA,'star',sym_color='orange',sym_size=2.5,/sym_filled) ;p = plot(lon_ind,lat_ind,symbol='o',/sym_filled,sym_size=0.25,linestyle=6,color='blue',/overplot) ;t = text(213,-6,'NOAA',color='red',font_size=10,/data) ;t = text(213,-12,'AGAGE',color='orange',font_size=10,/data) ;t = text(295,33,'SEAC$^4$RS ER-2',color='blue',font_size=10,/data) ;;t = text(250,55,'CarbonTracker',color='purple',font_size=10,/data) ;;t = text(250,51,'NAM',color='purple',font_size=10,/data) ;;t = text(250,-6,'CarbonTracker',color='green',font_size=10,/data) ;;t = text(250,-12,'Global',color='green',font_size=10,/data) ;map['Latitudes'].LABEL_ANGLE = 90 ;map['Longitudes'].LABEL_ANGLE = 0 ;map.save,dir+'Plots/SEAC4RS/Map_surface_sites.png' ;map.save,dir+'Plots/SEAC4RS/Map_surface_sites.pdf' end