; ; Calibrates, usually via regression, various NH and quasi-NH records ; against NH or quasi-NH seasonal or annual temperatures. ; ; Specify period over which to compute the regressions (stop in 1960 to avoid ; the decline that affects tree-ring density records) ; perst=1881. peren=1960. thalf=10. ; filter to use to give extra hi & lo pass info ; ; Select season of the temperature data against which to calibrate ; if n_elements(doseas) eq 0 then doseas=0 ; 0=annual, 1=Apr-Sep, 2=Oct-Mar seasname=['annual','aprsep','octmar'] ; ; Select record to calibrate ; if n_elements(dorec) eq 0 then dorec=0 recname=['esper','jones','mann','3tree','briffa','iwarm','mj2000','crowley03','rutherford04'] print,recname[dorec] doland=0 ; for Mann et al. only, 0=use NH mean, 1=use land>20N mean ; if recname[dorec] eq 'mann' then begin if doland eq 0 then openw,1,'recon_mannNH.out'+seasname[doseas] $ else openw,1,'recon_mannLN20.out'+seasname[doseas] endif else begin openw,1,'recon_'+recname[dorec]+'.out'+seasname[doseas] endelse ; multi_plot,nrow=2 if !d.name eq 'X' then window,ysize=800 ; ; Compute the >20N land instrumental temperature timseries and filter ; datst=1860 daten=1980 print,'Reading temperatures' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2_18512001.mon.nc') tsmon=crunc_rddata(ncid,tst=[datst,0],tend=[daten,11],grid=gtemp) ncdf_close,ncid nmon=12 ntemp=gtemp.nt nyrtemp=ntemp/nmon yrtemp=reform(gtemp.year,nmon,nyrtemp) yrtemp=reform(yrtemp(0,*)) ; ; Compute the northern hemisphere >20N land series ; ; First extract >20N rows kl=where(gtemp.y gt 20.) ylat=gtemp.y(kl) tsnorth=tsmon(*,kl,*) ; Compute latitude-weighted mean nhmon=globalmean(tsnorth,ylat) ; Compute seasonal/annual mean nhmon=reform(nhmon,nmon,nyrtemp) case doseas of 0: lvy=mkseason(nhmon,0,11,datathresh=6) ; could try 9,8 (Oct-Sep annual)! 1: lvy=mkseason(nhmon,3,8,datathresh=3) 2: lvy=mkseason(nhmon,9,2,datathresh=3) endcase ; ; Filter it ; filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan ; ; Read in record and filter ; case recname[dorec] of 'esper': begin openr,2,'/cru/u2/f055/data/paleo/esper2002/esper.txt' readf,2,nyr headst='' readf,2,headst rawdat=fltarr(7,nyr) readf,2,rawdat close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(1,*)) end 'jones': begin openr,2,'../tree5/phil_nhrecon.dat' nyr=992 rawdat=fltarr(4,nyr) readf,2,rawdat,format='(I5,F7.2,I3,F7.2)' close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(3,*)) end 'mann': begin if doland eq 0 then begin openr,2,'../tree5/mann_nhrecon1000.dat' nyr=981 rawdat=fltarr(2,nyr) readf,2,rawdat ;,format='(I6,F11.7)' close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(1,*)) endif else begin openr,2,'../tree5/mannarea_all.dat' nyr=981 rawdat=fltarr(11,nyr) headdat=' ' readf,2,headdat readf,2,rawdat ;,format='(I6,F11.7)' close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(10,*)) endelse end '3tree': begin openr,2,'../tree6/tornyamataim.ave' readf,2,nnn rawdat=fltarr(2,nnn) readf,2,rawdat close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(1,*)) end 'briffa': begin restore,filename='/cru/u2/f055/tree6/bandtempNH_calmultipcr.idlsave' x=yrmxd densall=prednh end 'iwarm': begin ; use warm-season instrumental series as the predictor! x=yrtemp densall=mkseason(nhmon,3,8,datathresh=3) end 'mj2000': begin openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/mann03_orig.dat' readf,2,ncol readf,2,icol readf,2,nyr readf,2,nhead headst=strarr(nhead) rawdat=fltarr(ncol,nyr) readf,2,headst readf,2,rawdat close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(icol,*)) end 'crowley03': begin openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/crowley03_orig.dat' readf,2,ncol readf,2,icol readf,2,nyr readf,2,nhead headst=strarr(nhead) rawdat=fltarr(ncol,nyr) readf,2,headst readf,2,rawdat close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(icol,*)) end 'rutherford04': begin openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/rutherford04_orig.dat' readf,2,ncol readf,2,icol readf,2,nyr readf,2,nhead headst=strarr(nhead) rawdat=fltarr(ncol,nyr) readf,2,headst readf,2,rawdat close,2 x=reform(rawdat(0,*)) densall=reform(rawdat(icol,*)) end endcase ; oopsx=x oopsy=densall ; kl=where((x ge datst) and (x le daten)) x=x(kl) densall=densall(kl) if total(abs(yrtemp-x)) ne 0 then message,'Incompatible years' y1=densall filter_cru,thalf,tsin=y1,tslow=ylow1,tshigh=yhi1,/nan ; ; Now correlate and regress them ; printf,1,'Correlations and regression coefficients for '+seasname[doseas] keeplist=where(finite(y1+lvy) and (x ge perst) and (x le peren),nkeep) x=x(keeplist) ts1=y1(keeplist) ts2=lvy(keeplist) r1=correlate(ts1,ts2) dum=linfit(ts1,ts2,sigma=err_sigma) c1=dum printf,1,'Full timeseries:',r1,c1 plot,x,ts2 oplot,x,ts1*c1[1]+c1[0],thick=2 ; tsa=ts1(0:nkeep-2) tsb=ts1(1:nkeep-1) ;mxd_ar1=correlate(tsa,tsb) mxd_ar1=a_correlate(ts1,[-1]) tsa=ts2(0:nkeep-2) tsb=ts2(1:nkeep-1) ;nht_ar1=correlate(tsa,tsb) nht_ar1=a_correlate(ts2,[-1]) printf,1,'AR1 for MXD and NHEMI:',mxd_ar1,nht_ar1 ; tspred=c1(0)+c1(1)*ts1 tswant=ts2 plot,tspred,ts2,psym=1,$ xtitle='Scaled Esper et al. anomaly (!Uo!NC)',$ ytitle='Northern Hemisphere temperature anomaly (!Uo!NC)',$ /xstyle,xrange=[-0.6,0.3],$ /ystyle,yrange=[-0.6,0.3] oplot,!x.crange,[0.,0.],linestyle=1 oplot,[0.,0.],!y.crange,linestyle=1 dum=linfit(tspred,ts2,sigma=se) oplot,!x.crange,!x.crange*dum(1)+dum(0) oplot,!x.crange,!x.crange*dum(1)+dum(0)+se(0) oplot,!x.crange,!x.crange*dum(1)+dum(0)-se(0) oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0) oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0) oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)+se(0) oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)-se(0) oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)-se(0) oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)+se(0) ; ts1=yhi1(keeplist) ts2=lvyhi(keeplist) r2=correlate(ts1,ts2) dum=linfit(ts1,ts2) c2=dum printf,1,'High-pass :',r2,c2 ; ts1=ylow1(keeplist) ts2=lvylow(keeplist) r3=correlate(ts1,ts2) dum=linfit(ts1,ts2) c3=dum printf,1,'Low-pass :',r3,c3 ; ; Now compute the rms error between the reconstruction and the original ; timeseries ; tserr=tswant-tspred rmserr=sqrt( total( tserr^2 ) / float(n_elements(tserr)) ) printf,1,'RMS error between land>20Napr-sep temperature and Esper et al. reconstruction' printf,1,rmserr printf,1,'Uncertainties surrounding regression coefficients' printf,1,err_sigma ; printf,1,' ' printf,1,'Computations carried out over the period ',perst,peren printf,1,' ' printf,1,'To separate low and high frequency components, a gaussian weighted' printf,1,'filter was used with a half-width (years) of ',thalf ; close,1 ; end