; ; Computes regressions on full, high and low pass Esper et al. (2002) series, ; anomalies against full NH temperatures and other series. ; CALIBRATES IT AGAINST THE LAND-ONLY TEMPERATURES NORTH OF 20 N ; ; Specify period over which to compute the regressions (stop in 1960 to avoid ; the decline ; perst=1881. peren=1960. doseas=0 ; 0=annual, 1=Apr-Sep, 2=Oct-Mar seasname=['annual','aprsep','octmar'] ; openw,1,'recon_esper.out'+seasname[doseas] thalf=10. ; ; Compute the >20N land instrumental temperature timseries and filter ; print,'Reading temperatures' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_jones_18512000.mon.nc') tsmon=crunc_rddata(ncid,tst=[1860,0],tend=[1992,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) 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 Esper timeseries and filter ; 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,*)) kl=where((x ge 1860) and (x le 1992)) 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 ; 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 ; multi_plot,nrow=2 if !d.name eq 'X' then window,ysize=800 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