; ; Computes regressions on full, high and low pass MEAN timeseries of MXD ; anomalies against full NH temperatures. ; THIS IS FOR THE Mann et al. reconstruction ; CALIBRATES IT AGAINST THE LAND-ONLY TEMPERATURES NORTH OF 20 N ; IN FACT, I NOW HAVE AN ANNUAL LAND-ONLY NORTH OF 20N VERSION OF MANN, ; SO I CAN CALIBRATE THIS TOO - WHICH MEANS I'm ONLY ALTERING THE SEASON ; doland=1 ; 0=use Mann NH 1=use Mann land north of 20N ; ; Specify period over which to compute the regressions (stop in 1940 to avoid ; the decline ; perst=1881. peren=1960. ; openw,1,'recon'+string(doland+2,format='(I1)')+'_mann.out' 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_18511998.mon.nc') tsmon=crunc_rddata(ncid,tst=[1860,0],tend=[1990,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 Apr-Sep mean nhmon=reform(nhmon,nmon,nyrtemp) lvy=mkseason(nhmon,3,8,datathresh=3) ; ; Filter it ; filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan ; ; Read in Mann timeseries and filter ; 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 kl=where((x ge 1860) and (x le 1990)) 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 between normalised timeseries' 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 age-banded MXD 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 NH temperature and age-banded-MXD reconstruction' printf,1,rmserr printf,1,'Uncertainties surrounding regression coefficients' printf,1,err_sigma ; ; Now do some example reconstructions (also for 1912 & 1884 we know the ; instrumental values) ; printf,1,' ' ele1=where(x eq 1884) printf,1,'Instrumental NH temp for 1884 is',tswant(ele1(0)) ele1=where(x eq 1912) printf,1,'Instrumental NH temp for 1912 is',tswant(ele1(0)) printf,1,' ' ; exyr=[1601,1816,1641,1453,1912,1884] exval=[-6.9026,-4.3256,-4.3124,-4.2402,-3.3282,-2.8874] nex=n_elements(exyr) printf,1,'Reconstructions with error bounds' for i = 0 , nex-1 do begin xx=exval(i) y1=xx*(c1(1)+err_sigma(1))+c1(0)-err_sigma(0)-rmserr y2=xx*c1(1)+c1(0) y3=xx*(c1(1)-err_sigma(1))+c1(0)+err_sigma(0)+rmserr printf,1,exyr(i),y1,y2,y3 endfor ; 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