;
; Calibrates, usually via regression, various NH and quasi-NH records 
; against NH or quasi-NH seasonal or annual temperatures.
;
; Specify period over which to compute the regressions (stop in 1960 to avoid
; the decline that affects tree-ring density records)
;
perst=1881.
peren=1960.
thalf=10.     ; filter to use to give extra hi & lo pass info
;
; 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']
;
; Select record to calibrate
;
if n_elements(dorec) eq 0 then dorec=0
recname=['esper','jones','mann','3tree','briffa','iwarm','mj2000','crowley03','rutherford04']
print,recname[dorec]
doland=0      ; for Mann et al. only, 0=use NH mean, 1=use land>20N mean
;
if recname[dorec] eq 'mann' then begin
  if doland eq 0 then openw,1,'recon_mannNH.out'+seasname[doseas] $
                 else openw,1,'recon_mannLN20.out'+seasname[doseas]
endif else begin
  openw,1,'recon_'+recname[dorec]+'.out'+seasname[doseas]
endelse
;
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=800
;
; Compute the >20N land instrumental temperature timseries and filter
;
datst=1860
daten=1980
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
;
; Filter it
;
filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan
;
; Read in record and filter
;
case recname[dorec] of
  'esper': begin
    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,*))
  end
  'jones': begin
    openr,2,'../tree5/phil_nhrecon.dat'
    nyr=992
    rawdat=fltarr(4,nyr)
    readf,2,rawdat,format='(I5,F7.2,I3,F7.2)'
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(3,*))
  end
  'mann': begin
    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
  end
  '3tree': begin
    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,*))
  end
  'briffa': begin
    restore,filename='/cru/u2/f055/tree6/bandtempNH_calmultipcr.idlsave'
    x=yrmxd
    densall=prednh
  end
  'iwarm': begin     ; use warm-season instrumental series as the predictor!
    x=yrtemp
    densall=mkseason(nhmon,3,8,datathresh=3)
  end
  'mj2000': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/mann03_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
  'crowley03': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/crowley03_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
  'rutherford04': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/rutherford04_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
endcase
;
oopsx=x
oopsy=densall
;
kl=where((x ge datst) and (x le daten))
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
plot,x,ts2
oplot,x,ts1*c1[1]+c1[0],thick=2
;
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
;
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