pro mp ; All the data used is from the WDC for paleoclimatology, http://www.ncdc.noaa.gov/paleo/icecore.html ;CO2: WDC PALEO CONTRIBUTION SERIES CITATION: ; Lüthi, D., et al.. 2008. ; EPICA Dome C Ice Core 800KYr Carbon Dioxide Data. ; IGBP PAGES/World Data Center for Paleoclimatology ; Data Contribution Series # 2008-055. ; NOAA/NCDC Paleoclimatology Program, Boulder CO, USA. ; IGBP PAGES/WDCA CONTRIBUTION SERIES NUMBER: 2008-055 ;CH4: WDC PALEO CONTRIBUTION SERIES CITATION: ; Loulergue, L., et al.. 2008. ; EPICA Dome C Ice Core 800KYr Methane Data. ; IGBP PAGES/World Data Center for Paleoclimatology ; Data Contribution Series # 2008-054. ; NOAA/NCDC Paleoclimatology Program, Boulder CO, USA. ; IGBP PAGES/WDCA CONTRIBUTION SERIES NUMBER: 2008-054 ;N2O: WDC PALEO CONTRIBUTION SERIES CITATION: ; Schilt, A., et al. 2010. ; EPICA Dome C Ice Core 800KYr N2O Data. ; IGBP PAGES/World Data Center for Paleoclimatology ; Data Contribution Series # 2010-005. ; NOAA/NCDC Paleoclimatology Program, Boulder CO, USA. ; IGBP PAGES/WDCA CONTRIBUTION SERIES NUMBER: 2010-005 ;All datasets are combined records on the EDC3 timescale ; the three filenames for the three tracegases ngas=3 fname=strarr(ngas) fname=['edc-co2-2008.txt','edc-ch4-2008.txt','edc-n2o-2010-800k.txt'] ; size of the header for each file nheader=intarr(ngas) nheader=[770,154,138] header='' ; sdata is the dimensions of the data to be read in for each gas sdata=intarr(ngas,2) sdata(0,*)=[2,1096] sdata(1,*)=[4,2103] sdata(2,*)=[4,911] ; all data is the longest record alldata=max(sdata(*,1)) ; ndata is the column index for the year and gas data ndata=intarr(ngas,2) ndata(0,*)=[0,1] ndata(1,*)=[1,2] ndata(2,*)=[2,3] ; plotrange is the y-axis range plotrange=fltarr(2,ngas) plotrange(*,0)=[240,290] plotrange(*,1)=[400,900] plotrange(*,2)=[200,300] ; namegas is the name for the plot title namegas=strarr(ngas) namegas=['Carbon Dioxide CO2','Methane CH4','Nitrous Oxide N2O'] gnamef=strarr(ngas) gnamef=['co2','ch4','n2o'] ; datatime=fltarr(alldata,ngas) datagas=fltarr(alldata,ngas) ; xyranges for the two plots xxrange=fltarr(2,2) xxrange(*,0)=[0,11] xxrange(*,1)=[114,133] ; timeslice names tnamef=strarr(2) tnamef=['hol','lig'] ; nt is the time resolution of the output pmip data. 0.1 = 100 years ; times is the array of times nt=0.1 nstart=0 nfinish=133 times=nstart+nt*findgen( 1+(nfinish-nstart)/nt) ; the pmip transients tpmip=fltarr(2,2) tpmip(*,0)=[0,10] tpmip(*,1)=[115,132] ; the pmip time slices ntime=4 timeslice=fltarr(ntime) timeslice=[115,125,128,130] ; loop over each gas/file: for g=0,ngas-1 do begin ; open the file for reading close,1 openr,1,fname(g) for h=0,nheader(g)-1 do begin readf,1,header endfor lline=fltarr(sdata(g,0)) ; read in a line at a time. Convert years to kyr. for l=0,sdata(g,1)-1 do begin readf,1,lline datatime(l,g)=lline(ndata(g,0))/1000.0 datagas(l,g)=lline(ndata(g,1)) endfor endfor close,1 ninterp=3 aa=size(times) co2interp=fltarr(aa(1),ngas,ninterp) loadct,39 for g=0,ngas-1 do begin ; OK, now do the interpolation with standard linear interpolation ; 0 is linear interpolation direct from data ; 1 applies a smoothing (10-point window) ; 2 is constant co2interp(*,g,0)=interpol([datagas(0,g),datagas(0:sdata(g,1)-1,g)],[0,datatime(0:sdata(g,1)-1,g)],times) co2interp(*,g,1)=smooth(co2interp(*,g,0),10,/edge_truncate,/nan) co2interp(*,0,2)=280.0 co2interp(*,1,2)=650.0 co2interp(*,2,2)=270.0 ; and finally plot it... set_plot,'ps' for t=0,1 do begin !P.FONT=0 device,filename='pmipgas_'+gnamef(g)+'_'+tnamef(t)+'_new.eps',/color,/encapsulate,xsize=25,ysize=15,set_font='Helvetica' plot,datatime(0:sdata(g,1)-1,g),datagas(0:sdata(g,1)-1,g),xrange=xxrange(*,t),yrange=plotrange(*,g),xstyle=1,ystyle=1,/nodata,title=namegas(g) oplot,datatime(0:sdata(g,1)-1,g),datagas(0:sdata(g,1)-1,g),psym=7,symsize=0.5 for ni=0,ninterp-2 do begin oplot,times(where(times eq tpmip(0,t)):where(times eq tpmip(1,t))),co2interp(where(times eq tpmip(0,t)):where(times eq tpmip(1,t)),g,ni),color=80+(ni*80),thick=5 for tt=0,ntime-1 do begin if (where(times eq timeslice(tt)) ne -1) then begin oplot,times(where(times eq timeslice(tt))),co2interp(where(times eq timeslice(tt)),g,ni),color=80+(ni*80),psym=6,symsize=3 endif endfor ; end tt endfor ; end ni device,/close endfor ; end t endfor ; end g close,2 openw,2,'pmip_hol_lig_gases.txt' ; for each pmip time slice, which interpolation routine do we want? ; use direct for LIG my_interp_pmip=intarr(ntime) my_interp_pmip=[0,0,0,0] ; for each pmip transient, which interpolation routine do we want? ; use direct for LIG ; use smoothesd for holocene my_interp_pmipt=intarr(1) my_interp_pmipt=[1,0] my_pmipname=strarr(2) my_pmipname=['Holocene','LIG'] printf,2,'' printf,2,'CO2 CH4 and N2O for PMIP snapshots:' printf,2,'time (kyr BP)','co2 (ppm) ','ch4 (ppb) ','n2o (ppb) ' for t=0,ntime-1 do begin co2=fix(co2interp(where(times eq timeslice(t)),0,my_interp_pmip(t))+0.5) ch4=fix(co2interp(where(times eq timeslice(t)),1,my_interp_pmip(t))+0.5) n2o=fix(co2interp(where(times eq timeslice(t)),2,my_interp_pmip(t))+0.5) printf,2,strtrim(timeslice(t),2)+' '+strtrim(co2,2)+' '+strtrim(ch4,2)+' '+strtrim(n2o,2) endfor printf,2,'' printf,2,'CO2 CH4 and N2O for PMIP transients:' printf,2,'time (kyr BP)','co2 (ppm) ','ch4 (ppb) ','n2o (ppb) ' for t=0,1 do begin printf,2,'' printf,2,my_pmipname(t) ni=my_interp_pmipt(t) s=size(times(where(times eq tpmip(0,t)):where(times eq tpmip(1,t)))) a=times(where(times eq tpmip(0,t)):where(times eq tpmip(1,t))) b=co2interp(where(times eq tpmip(0,t)):where(times eq tpmip(1,t)),*,ni) for x=0,s(1)-1 do begin printf,2,a(x),b(x,0),b(x,1),b(x,2) endfor endfor close,2 stop end