pro create_tracer_entry_time_series

dir = '/Volumes/Data_Gator/Data/Surface_Trace_Gas/'

;  Read in the tropical N2O time series.
restore,dir+'N2O_tropical_time_series.sav'
nt_n2o = n_elements(tropical_n2o_time)
;tropical_n2o_time += 2./12.

;  Read in the tropical CO2 time series.
restore,dir+'CO2_tropical_time_series.sav'
nt_co2 = n_elements(tropical_co2_time)
;tropical_co2_time += 2./12.

;  Read in the tropical CH4 time series.
restore,dir+'CH4_tropical_time_series.sav'

;  Read in the tropical SF6 time series.
restore,dir+'SF6_tropical_time_series.sav'
nt_sf6 = n_elements(sf6_time)
dt = sf6_time[1]-sf6_time[0]
;sf6_time += 2./12.

time_entry = tropical_n2o_time
nt = n_elements(time_entry)
print,tropical_n2o_time[0],tropical_n2o_time[-1],tropical_co2_time[0],tropical_co2_time[-1],tropical_ch4_time[0],tropical_ch4_time[-1],sf6_time[0],sf6_time[-1]
ti_co2 = where(tropical_co2_time ge time_entry[0] and tropical_co2_time le time_entry[-1])
ti_sf6 = where(sf6_time ge time_entry[0]-1e-2 and sf6_time le time_entry[-1]+1e-2)


;  Set up the idealized Green's function age spectra for a full range of mean ages in both the tropics and midlats.
dt_spectra = 1./12.
yrmax_spectra = 25.
nt_spectra = yrmax_spectra/dt_spectra+1
age_spectra_years = findgen(nt_spectra)*dt_spectra + 0.0005
age_spectra_years[0] = 0.001

;  Stratospheric ratios based on Fritsch et al., 2020.  Tropospheric ratios based on Boenisch et al., 2010.
ratios_trop = [findgen(6)*0.1+0.1,0.8,1,1.2,1.5,1.75,2,2.5]
ratios = [0.8,1.0,1.25]
nr = n_elements(ratios)
nr_trop = n_elements(ratios_trop)

;  Compute the Green's functions in the tropics and midlats for the full range of possible mean ages.
d_age = 0.1
max_age = 15.0
max_age_trop = 9.0
n_ages = max_age/d_age
n_ages_trop = max_age_trop/d_age
grn_mean_ages = findgen(n_ages)*d_age
grn_mean_ages[0] = 0.002
grn_mean_ages_trop = findgen(n_ages_trop)*d_age
grn_mean_ages_trop[0] = 0.002
grn = fltarr(nt_spectra,n_ages,nr) & delta_grn = fltarr(n_ages,nr) & grn_trop = fltarr(nt_spectra,n_ages_trop,nr_trop) & delta_grn_trop = fltarr(n_ages_trop,nr_trop)
for r = 0, nr-1 do for t = 0, n_ages-1 do begin
  grn[*,t,r] = 1./(2.*sqrt(ratios[r]*grn_mean_ages[t])*sqrt(!pi*(age_spectra_years/grn_mean_ages[t])^3)) * $
    exp((-(grn_mean_ages[t]^2*(age_spectra_years/grn_mean_ages[t]-1)^2)) / (4.*sqrt(ratios[r]*grn_mean_ages[t])^2*age_spectra_years/grn_mean_ages[t]))
  grn[*,t,r] /= total(grn[*,t,r])
  delta_grn[t] = sqrt(0.5 * total((age_spectra_years-grn_mean_ages[t])^2 * grn[*,t,r]))
endfor
for r = 0, nr_trop-1 do for t = 0, n_ages_trop-1 do begin
  grn_trop[*,t,r] = 1./(2.*sqrt(ratios_trop[r]*grn_mean_ages_trop[t])*sqrt(!pi*(age_spectra_years/grn_mean_ages_trop[t])^3)) * $
    exp((-(grn_mean_ages_trop[t]^2*(age_spectra_years/grn_mean_ages_trop[t]-1)^2)) / (4.*sqrt(ratios_trop[r]*grn_mean_ages_trop[t])^2*age_spectra_years/grn_mean_ages_trop[t]))
  grn_trop[*,t,r] /= total(grn_trop[*,t,r])
  delta_grn_trop[t] = sqrt(0.5 * total((age_spectra_years-grn_mean_ages_trop[t])^2 * grn_trop[*,t,r]))
endfor

lat_grid_entry = lat_grid_n2o

;  Convolve the Green's Fcns with the tpause time series.
ny = n_elements(lat_grid_entry)
aa = fltarr(n_ages,ny,nr,nt) & n2o_entry = aa & co2_entry = aa & ch4_entry = aa & sf6_entry = aa & n2o_entry_lat = aa & co2_entry_lat = aa & ch4_entry_lat = aa & sf6_entry_lat = aa
bb = fltarr(n_ages_trop,ny,nr_trop,nt) & n2o_entry_trop = bb & co2_entry_trop = bb & ch4_entry_trop = bb & sf6_entry_trop = bb & lat_entry_trop = fltarr(n_ages_trop,ny,nr_trop)
cc = fltarr(n_ages_trop,nr_trop,nt) & n2o_tropical_entry = cc & co2_tropical_entry = cc & ch4_tropical_entry = cc & sf6_tropical_entry = cc & frac_extratropics_trop = fltarr(n_ages_trop,nr_trop)

;  Time weighting of age spectra from equator and non-equator.
ti1 = 0.4
ti2 = 0.6
t1 = where(age_spectra_years le ti1,nt1)
t2 = where(age_spectra_years gt ti1 and age_spectra_years lt ti2,nt2)
t3 = where(age_spectra_years ge ti2,nt3)
linear1 = 1. - findgen(nt2)/nt2
gauss1 = exp(-(age_spectra_years[t2]-ti1)^2. / (ti2/30.*(ti2-ti1)*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
t_init = 26/dt_spectra

for t = t_init, nt-1 do for r = 0, nr-1 do for y = 0, ny-1 do for tt = 0, n_ages-1 do begin
  n2o_entry[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn[*,tt,r])*tropical_n2o_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn[*,tt,r])*tropical_n2o[lindgen(nt_spectra)+t-25./dt_spectra])
  co2_entry[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn[*,tt,r])*tropical_co2_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn[*,tt,r])*tropical_co2[lindgen(nt_spectra)+t-25./dt_spectra])
  ch4_entry[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn[*,tt,r])*tropical_ch4_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn[*,tt,r])*tropical_ch4[lindgen(nt_spectra)+t-25./dt_spectra])
  sf6_entry[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn[*,tt,r])*tropical_sf6_lat[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn[*,tt,r])*tropical_sf6[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra])
  n2o_entry_lat[tt,y,r,t] = total(reverse(grn[*,tt,r])*tropical_n2o_lat[lindgen(nt_spectra)+t-25./dt_spectra,y])
  co2_entry_lat[tt,y,r,t] = total(reverse(grn[*,tt,r])*tropical_co2_lat[lindgen(nt_spectra)+t-25./dt_spectra,y])
  ch4_entry_lat[tt,y,r,t] = total(reverse(grn[*,tt,r])*tropical_ch4_lat[lindgen(nt_spectra)+t-25./dt_spectra,y])
  sf6_entry_lat[tt,y,r,t] = total(reverse(grn[*,tt,r])*tropical_sf6_lat[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra,y])
endfor
chk = where(~finite(n2o_entry),nchk)
print,nchk
chk = where(~finite(co2_entry),nchk)
print,nchk
for t = t_init, nt-1 do for r = 0, nr_trop-1 do for y = 0, ny-1 do for tt = 0, n_ages_trop-1 do begin
  n2o_entry_trop[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])*tropical_n2o_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn_trop[*,tt,r])*tropical_n2o[lindgen(nt_spectra)+t-25./dt_spectra])
  co2_entry_trop[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])*tropical_co2_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn_trop[*,tt,r])*tropical_co2[lindgen(nt_spectra)+t-25./dt_spectra])
  ch4_entry_trop[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])*tropical_ch4_lat[lindgen(nt_spectra)+t-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn_trop[*,tt,r])*tropical_ch4[lindgen(nt_spectra)+t-25./dt_spectra])
  sf6_entry_trop[tt,y,r,t] = total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])*tropical_sf6_lat[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra,y] + $
    reverse(weight_tr)*reverse(grn_trop[*,tt,r])*tropical_sf6[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra])
endfor
for r = 0, nr_trop-1 do for y = 0, ny-1 do for tt = 0, n_ages_trop-1 do lat_entry_trop[tt,y,r] = total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])*replicate(lat_grid_entry[y],nt_spectra))
for r = 0, nr_trop-1 do for tt = 0, n_ages_trop-1 do frac_extratropics_trop[tt,r] = $
  total(reverse(weight_nh)*reverse(grn_trop[*,tt,r])) / total(reverse(weight_tr)*reverse(grn_trop[*,tt,r]) + reverse(weight_nh)*reverse(grn_trop[*,tt,r]))


;  Tropical entry only.
for t = t_init, nt-1 do for r = 0, nr_trop-1 do for tt = 0, n_ages_trop-1 do begin
  n2o_tropical_entry[tt,r,t] = total(reverse(grn_trop[*,tt,r])*tropical_n2o[lindgen(nt_spectra)+t-25./dt_spectra])
  co2_tropical_entry[tt,r,t] = total(reverse(grn_trop[*,tt,r])*tropical_co2[lindgen(nt_spectra)+t-25./dt_spectra])
  ch4_tropical_entry[tt,r,t] = total(reverse(grn_trop[*,tt,r])*tropical_ch4[lindgen(nt_spectra)+t-25./dt_spectra])
  sf6_tropical_entry[tt,r,t] = total(reverse(grn_trop[*,tt,r])*tropical_sf6[lindgen(nt_spectra)+t+ti_sf6[0]-25./dt_spectra])
endfor

save,time_entry,lat_grid_entry,grn_mean_ages,grn_mean_ages_trop,ratios,ratios_trop,n2o_entry,co2_entry,ch4_entry,sf6_entry,n2o_entry_lat,co2_entry_lat,ch4_entry_lat,sf6_entry_lat, $
  n2o_entry_trop,co2_entry_trop,ch4_entry_trop,sf6_entry_trop,grn,grn_trop,age_spectra_years,dt_spectra,yrmax_spectra,n2o_tropical_entry,co2_tropical_entry,ch4_tropical_entry, $
  sf6_tropical_entry,lat_entry_trop,frac_extratropics_trop,filename=dir+'Tracer_entry_time_series.sav'


ch4_surface_lat_norm = tropical_ch4_lat & n2o_surface_lat_norm = tropical_n2o_lat & co2_surface_lat_norm = tropical_co2_lat
for t = 0, nt-1 do begin
  ch4_surface_lat_norm[t,*] = reform(tropical_ch4_lat[t,*]) / tropical_ch4[t]
  n2o_surface_lat_norm[t,*] = reform(tropical_n2o_lat[t,*]) / tropical_n2o[t]
  co2_surface_lat_norm[t,*] = reform(tropical_co2_lat[t,*]) / tropical_co2[t]
endfor


colors = ['Blue','dark blue','Purple','Medium Purple','Maroon','Lavender','Magenta','Lime Green','Green','Orange','Red','Brown','Papaya Whip','Silver','Grey','Dark Grey']

t = 2
;p = plot(indgen(2),/nodata,xrange=[0,5])
;for r = 0, nr-1 do p = plot(age_spectra_years,grn[*,t,r],color=colors[r],thick=2,/overplot)
;;p = plot(age_spectra_years,grn_s[*,t],color='sky blue',thick=2,/overplot)
;;p = plot(age_spectra_years,grn_m[*,t],color='lime green',thick=2,/overplot)
;;p = plot(age_spectra_years,grn_trop2[*,t],color='pink',thick=2,/overplot)
;for r = 0, nr-1 do p = plot(age_spectra_years,weight_tr*grn[*,t,r],color=colors[r],/overplot)
;p = plot(age_spectra_years,weight_tr*grn_s[*,t],color='red',/overplot)
;p = plot(age_spectra_years,weight_tr*grn_m[*,t],color='green',/overplot)
;p = plot(age_spectra_years,weight_tr*grn_trop[*,t],color='purple',/overplot)
;p = plot(age_spectra_years,weight_tr*grn_trop2[*,t],color='grey',/overplot)

p = plot(indgen(2),/nodata,xrange=[0,6])
for r = 0, nr_trop-1 do p = plot(age_spectra_years,grn_trop[*,10,r],color=colors[r],thick=2,/overplot)

;p = plot(indgen(2),/nodata,xrange=[1.07,0.95],xtitle='Normalized CH$_4$',ytitle='Normalized N$_2$O',yrange=[1.005,0.995],xticklen=0.02,yticklen=0.02,title='CH$_4$ vs N$_2$O', $
;  font_size=11,dimensions=[700,400])
;for t = 660, 671 do p = plot(ch4_surface_lat_norm[t,*],n2o_surface_lat_norm[t,*],symbol='o',sym_size=1,/sym_filled,linestyle=6,color=colors[0:ny-1],/overplot)

;p = plot(indgen(2),/nodata,xrange=[0,1])
;p = plot(age_spectra_years,weight_tr,color='orange',/overplot)
;p = plot(age_spectra_years,weight_nh,color='red',/overplot)

end