pro make_tropical_ch4_time_series

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

;  Read in GMD marine boundary layer data for CH4.
dat = read_ascii(dir+'CH4/CH4_mbl_17S-17N.txt',data_start=72)
ch4_time = reform(dat.field1[3,*])
ch4_mbl = reform(dat.field1[4,*])
nt = n_elements(ch4_time)

;  Read in the CMIP6 surface time series.
dat = read_ascii(dir+'CMIP6_surface_tracer_time_series.txt',data_start=22)
cmip_time = reform(dat.field01[0,*]) + 0.5
cmip_ch4 = reform(dat.field01[2,*])

;  Read in latitude averaged mbl data from https://gml.noaa.gov/ccgg/obspack/our_products.php.
ny = 14
lat_grid_ch4 = indgen(ny)*10 - 60
ch4_mbl_tropics_lat = fltarr(nt,ny) & ratio = ch4_mbl_tropics_lat & ratio_seas = fltarr(12,ny)
for y = 0, ny-1 do begin
  if lat_grid_ch4[y] lt 0 then hem = 'S' else hem = 'N'
  lat_string = strcompress(string(abs(lat_grid_ch4[y])),/r)+hem
  dat = read_ascii(dir+'CH4/CH4_mbl_'+lat_string+'.txt',data_start=72)
  ch4_mbl_time = reform(dat.field1[3,*])
  ch4_mbl_tropics_lat[*,y] = reform(dat.field1[4,*])
  ratio[*,y] = ch4_mbl_tropics_lat[*,y] / ch4_mbl
  mon_fracs = ch4_mbl_time - fix(ch4_mbl_time)
  for t = 0., 11. do begin
    ti = where(mon_fracs ge t/12. and mon_fracs lt (t+1)/12.,nti)
    if nti gt 0 then ratio_seas[t,y] = mean(ratio[ti,y])
  endfor
endfor

;
;ii = interpol(indgen(nt_ch4),ch4_gv_time,ch4_mlo_time)
;ch4_gv_tropics_sm_interp = interpolate(ch4_gv_tropics_sm,ii)
;
;ti = where(ch4_mlo_time gt 2009)
;ch4_gv_tropics_sm_interp[ti] = 0.983 * ch4_mlo_sm[ti]
;
ch4_time_spo_core = findgen(44) + 1939.5
ch4_spo_core = [1090,1100,1100,1080,1085,1110,1115,1110,1110,1115,1130,1150,1130,1135,1160,1165,1170,1180,1190,1200,1220,1230,1240,1250,1260,1280,1290,1300,1315, $
  1340,1350,1370,1370,1410,1415,1420,1440,1450,1470,1500,1525,1550,1575,1600]

ch4_global_core = 1.02 * ch4_spo_core

ch4_time_etheridge = findgen(25)*2 + 1944.
ch4_etheridge = [1109,1120,1133,1147,1163,1182,1202,1224,1247,1272,1298,1326,1355,1386,1417,1449,1481,1514,1547,1580,1614,1644,1671,1694,1714]


cmip_diff = mean(ch4_mbl[where(ch4_time ge 1985 and ch4_time le 1995)]) - $
  mean(cmip_ch4[where(cmip_time ge 1985 and cmip_time le 1995)])

;ti = where(cmip_time lt ch4_time[0])
;ch4_mbl_time_o = [cmip_time[ti],ch4_time]
;ch4_mbl_o = [cmip_ch4[ti]+cmip_diff,ch4_mbl]

ti = where(ch4_time_etheridge lt ch4_time[0])
ch4_mbl_time_o = [ch4_time_etheridge[ti],ch4_time]
ch4_mbl_o = [ch4_etheridge[ti],ch4_mbl]

;ch4_mbl_time_o = [ch4_time_spo_core,ch4_time]
;ch4_mbl_o = [ch4_spo_core,ch4_mbl]

nt = round((2024.62 - 1945.) * 12. + 1)
dyear = 1./12
tropical_ch4_time = findgen(nt) * dyear + 1945 + 1./24
tropical_ch4_lat = fltarr(nt,ny)

tii = interpol(indgen(n_elements(ch4_mbl_time_o)),ch4_mbl_time_o,tropical_ch4_time)
tropical_ch4 = interpolate(ch4_mbl_o,tii)

tropical_ch4_sm = smooth(tropical_ch4,48,/edge_truncate)
deriv_yr = 12*sg_smooth(tropical_ch4_sm, NLEFT=12/2, NRIGHT=12/2-1, DERIV=1, DELTA=1.0, EDGE_PFCAST=1)
ch4_gr = 100 * deriv_yr / tropical_ch4_sm
ch4_gr_avg = mean(ch4_gr[-120:-24])

tii = where(tropical_ch4_time gt ch4_mbl_time_o[-1],nn)
for t = 0, nn-1 do tropical_ch4[tii[t]] = (1. + 1e-2*ch4_gr_avg) * tropical_ch4[tii[t]-12]

;  Extend latitudinal time series backward using seasonal cycle ratios.
ti = where(tropical_ch4_time lt ch4_mbl_time[0],nti)
ti2 = where(tropical_ch4_time ge ch4_mbl_time[0],nti2)
tii2 = interpol(indgen(n_elements(ch4_mbl_time)),ch4_mbl_time,tropical_ch4_time[ti2])
ti3 = where(tropical_ch4_time ge ch4_mbl_time_o[-1],nti3)
print,tropical_ch4_time[ti3]
for y = 0, ny-1 do begin
  for t = 0, nti-1 do tropical_ch4_lat[t,y] = ratio_seas[t mod 12,y] * tropical_ch4[t]
  tropical_ch4_lat[ti2,y] = interpolate(ch4_mbl_tropics_lat[*,y],tii2)
;  if nti3 gt 0 then for t = 0, nti3-1 do tropical_ch4_lat[ti3[t],y] = ratio_seas[t mod 1/dyear,y] * tropical_ch4[ti3[t]]
  if nti3 gt 0 then for t = 0, nti3-1 do tropical_ch4_lat[ti3[t],y] = (1. + 1e-2*ch4_gr_avg) * tropical_ch4_lat[ti3[t]-12,y]
endfor

;
save,tropical_ch4_time,tropical_ch4,tropical_ch4_lat,filename=dir+'CH4_tropical_time_series.sav'

gd = where(tropical_ch4_time ge 1945,ngd)

;  Write out CH4 in a text file.
;openw,1,dir+'CH4_tropical_surface_time_series.txt'
;printf,1,'Tropical CH4 mixing ratios'
;printf,1,'1945-1984: From Etheridge et al. (1998) based on scaled South Pole firn air.'
;printf,1,'1984-2023: MBL data set from NOAA (https://gml.noaa.gov/ccgg/mbl/).
;printf,1
;printf,1,'        Year       CH4 (ppb)'
;for t = 0, ngd-1 do printf,1,tropical_ch4_time[gd[t]],tropical_ch4[gd[t]],format='(2F13.3)'
;close,1


colors = ['sky blue','blue','purple','brown','pink','orange','red','lime green','green','tan','gold','medium orchid','lavender','grey']

;p = plot(co2_mbl_time,co2_mbl_tropics,thick=2,color='blue',xrange=[1940,2023],yrange=[315,420])

;p = plot(indgen(2),/nodata,xrange=[1980,2023])
;for y = 0, ny-1 do p = plot(ch4_mbl_time,ratio[*,y],thick=2,color=colors[y],/overplot)

p = plot(indgen(2),/nodata,xrange=[2020,2025],yrange=[1750,2050])
for y = 0, ny-1 do p = plot(tropical_ch4_time,tropical_ch4_lat[*,y],color=colors[y],/overplot)
p = plot(cmip_time,cmip_ch4,symbol='o',linestyle=6,color='red',/overplot)
p = plot(tropical_ch4_time,tropical_ch4,thick=2,color='magenta',/overplot)

;p = plot(indgen(2),/nodata,xrange=[2010,2024])
;p = plot(tropical_ch4_time,ch4_gr,thick=2,color='blue',/overplot)

;p = plot(tropical_ch4_time,tropical_ch4,thick=2,color='blue',xrange=[1960,2023],yrange=[1200,1950])

p = plot(indgen(2),/nodata,xrange=[1945,2025],yrange=[1000,2050])
;p = plot(ch4_gv_time,ch4_gv_tropics,thick=2,color='blue',/overplot)
;p = plot(ch4_gv_time,ch4_gv_tropics_sm,thick=2,color='green',/overplot)
;p = plot(ch4_mlo_time,ch4_gv_tropics_sm_interp,thick=2,color='lime green',/overplot)
;p = plot(ch4_mlo_time,ch4_mlo,thick=2,color='red',/overplot)
;p = plot(ch4_mlo_time,ch4_mlo_sm,thick=2,color='orange',/overplot)
;p = plot(ch4_spo_time,ch4_spo,thick=2,color='brown',/overplot)
;p = plot(ch4_spo_time,ch4_spo_sm,thick=2,color='khaki',/overplot)
p = plot(tropical_ch4_time,tropical_ch4,thick=2,color='magenta',/overplot)
p = plot(cmip_time,cmip_ch4,symbol='o',linestyle=6,color='red',/overplot)
;p = plot(ch4_time_spo_core,ch4_spo_core,thick=2,color='purple',/overplot)
p = plot(ch4_time_etheridge,ch4_etheridge,symbol='o',linestyle=6,color='sky blue',/overplot)
;p = plot(ch4_global_time,ch4_global,thick=2,color='plum',/overplot)

;p = plot(ch4_mlo_time,ch4_gv_tropics_sm_interp/ch4_mlo_sm,thick=2,color='green',yrange=[0.95,1.05])

;p = plot(ch4_mlo_time,ch4_gv_tropics_sm_interp/ch4_spo_sm,thick=2,color='green',yrange=[0.95,1.05])

;p = plot(indgen(2),/nodata)
;for y = 0, ny-1 do p = plot(indgen(12),ratio_seas[*,y],thick=2,color=colors[y],/overplot)

end