;
; Computes regressions on full, high and low pass MEAN timeseries of MXD
; anomalies against full NH temperatures.
; THIS IS FOR THE OVERPECK CIRCUM-ARCTIC RECORD
; CALIBRATES IT AGAINST THE LAND-ONLY TEMPERATURES NORTH OF 20 N
; *** Complicated because Overpeck record is five-yr-means only ***
;
; Specify period over which to compute the regressions (stop in 1940 to avoid
; the decline
;
perst=1880.
peren=1960.
;
; 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']
;
openw,1,'recon_overpeck.out'+seasname[doseas]
;
; Compute the >20N land instrumental temperature timseries and filter
;
datst=1868
daten=1992
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
;
; Compute five-yr means (1990 = 1988-1992 etc.)
;
declen=5
ndec=nyrtemp/declen
lvy=reform(lvy,declen,ndec)
lvy=total(lvy,1)/float(declen)
yrtemp=reform(yrtemp,declen,ndec)
yrtemp=reform(yrtemp(2,*))
;
; Read in Overpeck series
;
openr,2,'../tree6/overpeck_arctic.dat'
headst=strarr(2)
readf,2,headst
readf,2,nnn
rawdat=fltarr(2,nnn)
readf,2,rawdat
close,2
x=reform(rawdat(0,*))
densall=reform(rawdat(1,*))
kl=where((x ge 1870) and (x le 1990))
x=reverse(x(kl))
densall=reverse(densall(kl))
if total(abs(yrtemp-x)) ne 0 then message,'Incompatible years'
y1=densall
;
; 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)
;
printf,1,' '
printf,1,'Computations carried out over the period ',perst,peren
printf,1,'On 5-yr-means'
;
close,1
;
end