;
; Computes regressions on full, high and low pass MEAN timeseries of MXD
; anomalies against full NH temperatures.
; THIS IS FOR THE AVERAGE OF TORNE-YAMAL-TAIMYR
; 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_tyt.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 MXD timeseries and filter
;
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,*))
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