pro co2_surface_data

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

nco2_mlo = 190
co2_time = findgen(nco2_mlo)*1./12. + 2003.0 + 1./24.
nco2_ext = 430
co2_time_ext = findgen(nco2_ext)*1./12. + 1983.0 + 1./24.

co2_lats = [-13.75,19.5,25.7,40.0,46.0,71.0]
ny_co2 = n_elements(co2_lats)

co2_surface_flask = fltarr(nco2_ext,ny_co2)

data = read_ascii(dir+'CO2_lef.txt',data_start=69)
lef_yr = reform(data.field01[1,*])
lef_mon = reform(data.field01[2,*])
lef_co2_raw = reform(data.field01[11,*])
lef_time_raw = lef_yr + (lef_mon-1.)/12. + 1./24.
lef_time = lef_time_raw[uniq(lef_time_raw)]
chk = where(lef_time_raw ge 2013 and lef_time_raw lt 2014)
print,lef_time_raw[chk],lef_mon[chk],lef_co2_raw[chk]
nlef = n_elements(lef_time)
co2_lef_mon = replicate(!values.f_nan,nlef)
chk = where(lef_time ge 2013 and lef_time lt 2014)
for t = 0, nlef-1 do begin
  ti = where(lef_time_raw eq lef_time[t] and lef_co2_raw gt 0 and lef_co2_raw lt 450,nti)
  if nti gt 0 then co2_lef_mon[t] = mean(lef_co2_raw[ti])
endfor
print,lef_time[chk],co2_lef_mon[chk]
ti = where(lef_time ge co2_time[0]-0.01 and lef_time le co2_time[-1]+0.01)
co2_lef = co2_lef_mon[ti]
co2_surface_flask[240:-1,4] = co2_lef
for t = 239, 0, -1 do co2_surface_flask[t,4] = 0.9925 * co2_surface_flask[t+24,4]

data = read_ascii(dir+'CO2_sgp_monthly.txt',data_start=96)
sgp_yr = reform(data.field1[1,*])
sgp_mon = reform(data.field1[2,*])
co2_sgp_mon = reform(data.field1[3,*])
sgp_time = sgp_yr + (sgp_mon-1.)/12. + 1./24.
ti = where(sgp_time ge co2_time[0]-0.01 and sgp_time le co2_time[-1]+0.01)
co2_sgp = co2_sgp_mon[ti]
co2_surface_flask[240:-1,3] = co2_sgp
for t = 239, 0, -1 do co2_surface_flask[t,3] = 0.9925 * co2_surface_flask[t+24,3]

data = read_ascii(dir+'CO2_key_monthly.txt',data_start=69)
key_yr = reform(data.field1[1,*])
key_mon = reform(data.field1[2,*])
co2_key_mon = reform(data.field1[3,*])
key_time = key_yr + (key_mon-1.)/12. + 1./24.
ti = where(key_time ge co2_time[0]-0.01 and key_time le co2_time[-1]+0.01)
co2_key = co2_key_mon[ti]
co2_surface_flask[240:-1,2] = co2_key
for t = 239, 0, -1 do co2_surface_flask[t,2] = 0.9925 * co2_surface_flask[t+24,2]

data = read_ascii(dir+'CO2_mlo_monthly.txt',data_start=96)
mlo_yr = reform(data.field1[1,*])
mlo_mon = reform(data.field1[2,*])
co2_mlo_mon = reform(data.field1[3,*])
mlo_time = mlo_yr + (mlo_mon-1.)/12. + 1./24.
ti = where(mlo_time ge co2_time[0]-0.01 and mlo_time le co2_time[-1]+0.01)
co2_mlo = co2_mlo_mon[ti]
co2_surface_flask[240:-1,1] = co2_mlo
for t = 239, 0, -1 do co2_surface_flask[t,1] = 0.9925 * co2_surface_flask[t+24,1]

data = read_ascii(dir+'CO2_smo_monthly.txt',data_start=70)
smo_yr = reform(data.field1[1,*])
smo_mon = reform(data.field1[2,*])
co2_smo_mon = reform(data.field1[3,*])
smo_time = smo_yr + (smo_mon-1.)/12. + 1./24.
ti = where(smo_time ge co2_time[0]-0.01 and smo_time le co2_time[-1]+0.01)
co2_smo = co2_smo_mon[ti]
co2_surface_flask[240:-1,0] = co2_smo
for t = 239, 0, -1 do co2_surface_flask[t,0] = 0.9925 * co2_surface_flask[t+24,0]

data = read_ascii(dir+'CO2_brw_monthly.txt',data_start=70)
brw_yr = reform(data.field1[1,*])
brw_mon = reform(data.field1[2,*])
co2_brw_mon = reform(data.field1[3,*])
brw_time = brw_yr + (brw_mon-1.)/12. + 1./24.
ti = where(brw_time ge co2_time[0]-0.01 and brw_time le co2_time[-1]+0.01)
co2_brw = co2_brw_mon[ti]
co2_surface_flask[240:-1,5] = co2_brw
for t = 239, 0, -1 do co2_surface_flask[t,5] = 0.9925 * co2_surface_flask[t+24,5]


;  Read in Carbon Tracker data.
filenames = file_search(dir+'CarbonTracker/CT2019B.molefrac_glb*',count=nf)
filenames_nam = file_search(dir+'CarbonTracker/CT2019B.molefrac_nam*')
yrs = strmid(filenames,9,4,/reverse_offset)
mons = strmid(filenames,4,2,/reverse_offset)
yrmons = yrs+mons
ct_time = float(yrs) + float(mons-1.)/12. + 1./24.
ti = (where(co2_time_ext lt ct_time[0]))[-1]

ny_ct = 18
lat_ct = findgen(ny_ct)*10-85.

co2_ct_glb_zm = replicate(!values.f_nan,nco2_ext,ny_ct) & co2_ct_nam_zm = co2_ct_glb_zm & gph_nam_zm = co2_ct_glb_zm & blh_nam_zm = co2_ct_glb_zm & pres_nam_zm = co2_ct_glb_zm
for t = 0, nf-6 do begin
  ncdf_get,filenames[t],['longitude','latitude','pressure','time_components','decimal_date','blh','gph','co2'],ct,/quiet
  lon = ct['longitude','value']
  nx = (ct['longitude','dim_sizes'])[0]
  lat = ct['latitude','value']
  pres = ct['pressure','value']
  time = ct['time_components','value']
  blh = ct['blh','value']
  gph = ct['gph','value']
  co2_all = ct['co2','value']
  for y = 0, ny_ct-1 do begin
    yi = where(lat ge lat_ct[y]-5 and lat lt lat_ct[y]+5,nyi)
    co2_bl = fltarr(nx,nyi)
    for yy = 0, nyi-1 do for x = 0, nx-1 do begin
;      zi = where(gph[x,yi[yy],*]-gph[x,yi[yy],0] le blh[x,yi[yy]])
      zi = [0,1,2]
      co2_bl[x,yy] = mean(co2_all[x,yi[yy],zi])
    endfor
    co2_ct_glb_zm[ti+t+1,y] = mean(co2_bl)
  endfor

  ncdf_get,filenames_nam[t],['longitude','latitude','pressure','time_components','decimal_date','blh','gph','orography','co2'],ct,/quiet
  lon = ct['longitude','value']
  lat = ct['latitude','value']
  pres = 1e-2*ct['pressure','value']
  time = ct['time_components','value']
  blh = ct['blh','value']
  gph = ct['gph','value']
;  orography = ct['orography','value']
  co2_all = ct['co2','value']
  for y = 0, ny_ct-1 do begin
    yi = where(lat ge lat_ct[y]-5 and lat lt lat_ct[y]+5,nyi)
;    if t eq 200 then print,y,nyi
    if nyi gt 0 then begin
      xi = where(lon ge -112 and lon lt -85,nx)
      co2_bl = fltarr(nx,nyi) & gph_bl = co2_bl & blh_bl = co2_bl & pres_bl = co2_bl
      for yy = 0, nyi-1 do for x = 0, nx-1 do begin
;        zi = where(gph[xi[x],yi[yy],*]-gph[xi[x],yi[yy],0] le blh[xi[x],yi[yy]])
        zi = [0,1,2]
        co2_bl[x,yy] = mean(co2_all[xi[x],yi[yy],zi])
        gph_bl[x,yy] = gph[xi[x],yi[yy],3]-gph[xi[x],yi[yy],0]
        pres_bl[x,yy] = pres[xi[x],yi[yy],0]-pres[xi[x],yi[yy],-1]
        blh_bl[x,yy] = blh[xi[x],yi[yy]]
        if t eq 200 and y eq 12 then print,pres_bl[x,yy],gph_bl[x,yy],blh_bl[x,yy]
      endfor
      co2_ct_nam_zm[ti+t+1,y] = mean(co2_bl)
      gph_nam_zm[ti+t+1,y] = mean(gph_bl,/nan)
      blh_nam_zm[ti+t+1,y] = mean(blh_bl,/nan)
      pres_nam_zm[ti+t+1,y] = mean(pres_bl,/nan)
    endif
  endfor
endfor
print,pres_nam_zm[-30:-20,12],gph_nam_zm[-30:-20,12],blh_nam_zm[-30:-20,12]
for y = 0, ny_ct-1 do for t = ti, 0, -1 do begin
  co2_ct_glb_zm[t,y] = 0.9925 * co2_ct_glb_zm[t+24,y]
  co2_ct_nam_zm[t,y] = 0.9925 * co2_ct_nam_zm[t+24,y]
endfor

co2_ct_combo_zm = co2_ct_glb_zm
co2_ct_combo_zm[*,where(lat_ct ge 20 and lat_ct le 60)] = co2_ct_nam_zm[*,where(lat_ct ge 20 and lat_ct le 60)]

save,co2_time,co2_time_ext,nco2_mlo,nco2_ext,ny_co2,co2_lats,co2_surface_flask,lat_ct,ny_ct,co2_ct_glb_zm,co2_ct_nam_zm,co2_ct_combo_zm,filename=dir+'co2_surface_flask.sav'


colors = ['purple','blue','sky blue','lime green','green','orange','orange red','black','gray','sienna']

months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec','Jan']

p = plot(indgen(2),/nodata,xrange=[2013,2014],xtickname=months,xminor=0,xmajor=13,xtitle='2013',yrange=[383,408],ytitle='ppmv',title='GML Surface CO$_2$ Time Series',dimensions=[700,500], $
  margin=0.09,font_size=11)
p = plot([2013.67,2013.67],[383,408],linestyle=2,/overplot)
for y = 2, 5 do p = plot(co2_time_ext,co2_surface_flask[*,y],thick=3,color=colors[y+1],/overplot)
;for y = 9, 10 do p = plot(co2_time_ext,co2_ct_glb_zm[*,y],thick=3,color='medium orchid',/overplot)
;for y = 11, 14 do p = plot(co2_time_ext,co2_ct_nam_zm[*,y],thick=3,color=colors[y-11],/overplot)
for y = 11, 16 do p = plot(co2_time_ext,co2_ct_combo_zm[*,y],thick=3,color='medium orchid',/overplot)

;p = plot(indgen(2),/nodata,xrange=[1990,2014],xtitle='Year',yrange=[350,405],ytitle='ppmv',title='GML Surface CO$_2$ Time Series',dimensions=[700,500],margin=0.09,font_size=11)
;for y = 0, 1 do p = plot(co2_time_ext,co2_surface_flask[*,y],thick=3,color=colors[y+1],/overplot)
;for y = 7, 9 do p = plot(co2_time_ext,co2_ct_glb_zm[*,y],thick=3,color='medium orchid',/overplot)

end