; ; Computes regressions on full, high and low pass MEAN timeseries of MXD ; anomalies against full NH temperatures. ; THIS IS FOR THE AGE-BANDED (ALL BANDS) STUFF OF HARRY'S ; CALIBRATES IT AGAINST THE LAND-ONLY TEMPERATURES NORTH OF 20 N ; ; Specify period over which to compute the regressions (stop in 1940 to avoid ; the decline ; perst=1881. peren=1960. ; openw,1,'recon2_nhemi.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=[1997,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) ; ; Output it, then truncate at 1990 ; openw,5,'nhland20_amjjas.dat' printf,5,nyrtemp for iii = 0 , nyrtemp-1 do printf,5,yrtemp(iii),lvy(iii),format='(I4,F10.2)' close,5 kl=where(yrtemp le 1990,nyrtemp) yrtemp=yrtemp(kl) lvy=lvy(kl) ; ; Filter it ; filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan ; ; Read in MXD timeseries and filter ; restore,filename='densadj_all(330).idlsave' kl=where((x ge 1860) and (x le 1990)) x=x(kl) densall=densadj(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