;
; Plots long reconstructions (Age-banded, Jones, Mann(v2) and
; Yamal+Torne+Taimyr), and correlates them too!
;
doerr=0   ; 0=no errors, 1=Briffa age errors
;
; Have to restore the NEW age-banded results now to avoid overwriting key
; variables later
;
; For regional calibrated stuff
;restore,filename='bandtemp_calibrated.idlsave'
; Gets: calregts,nreg,nyr,regname,tempnyr,temptimey,tempregts,timey
;if regname(nreg-1) ne 'NH' then message,'Re-calibrate the age-banded!'
;newagetime=timey
;newagets=calregts(*,nreg-1)
;
; For allsites calibrated stuff
;restore,filename='bandalltemp_calibrated.idlsave'
; Gets: nyr,nhtit,timey,nhtemp
;if nhtit ne 'NH' then message,'Re-calibrate the age-banded!'
;newagetime=timey
;newagets=nhtemp
;
; For PCR calibrated stuff. The sm50 gets the smoothed series and appropriate
; erros - no need to smooth them!!!!  Except we want correlations between
; unsmoothed series, and between high-passed series, so we need to unsmoothed
; time series but the smoothed error bands!!
;
restore,filename='bandtempNH_calmultipcr.idlsave'
newagetime=yrmxd
newagets=prednh
restore,filename='bandtempNHsm50_calmultipcr.idlsave'
newagelow=prednh
newagese=predse
;restore,filename='bandtempNHsm50_calmultipcr_NSIBhug.idlsave'
; Gets: nyr,nhtit,yrmxd,prednh,fullnh,predse
if nhtit ne 'NH' then message,'Re-calibrate the age-banded!'
;
docorr=0        ; 1=compute correlations  0=don't
thalf=30.       ; filter for high or low-pass
;
lcol=[10,11,12,13,14,16,15]
lcol=replicate(15,7)
lthi=[4,2,4,9,9,4,6]
lsty=[0,0,2,9,9,1,0]
;
loadct,39
multi_plot,nrow=3,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=900
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,bold=0,font_size=16
endelse
def_1color,10,color='red'
def_1color,11,color='lpurple'
def_1color,12,color='mlblue'
def_1color,13,color='green'
def_1color,14,color='mdyellow'
def_1color,15,color='black'
def_1color,16,color='orange'
def_1color,18,color='white'
def_1color,22,color='vlblue'
def_smearcolor,fromto=[18,22]
;
ndo=7
nyr=2000
alltime=findgen(nyr)
alldo=fltarr(ndo,nyr)*!values.f_nan
namedo=['Jones ','Mann  ','Briffa','NChron','Overpeck','Crowley & Lowery',$
  'Instr']
iord=[0,1,2,5,6]
;
for jdo = 0 , n_elements(iord)-1 do begin
  ido=iord(jdo)
  fac=1.
  ;
  case ido of
    0: begin        ; Phil's reconstruction
      alltit="Jones' Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000
      peren=1992
      openr,1,'../tree5/phil_nhrecon.dat'
      nyr=992
      rawdat=fltarr(4,nyr)
      readf,1,rawdat,format='(I5,F7.2,I3,F7.2)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(3,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
      ; Convert from normalised values to deg C relative to 1961-1990
;      ts=ts*0.521 -0.1134       ; via Phil's variance matching (vs. NH all)
      ts=ts*0.3856-0.1112       ; via regression vs. NH land>20
    end
    1: begin        ; Mike Mann's reconstruction (now use his 1000 yr one)
      alltit="Mann's Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000
      peren=1980
      openr,1,'../tree5/mann_nhrecon1000.dat'
      nyr=981
      headdat=' '
      rawdat=fltarr(2,nyr)
;      readf,1,headdat
      readf,1,rawdat           ;,format='(I6,F11.7)'
      close,1
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
;      ts=ts(kl)-0.12            ; convert to oC wrt 1961-90 (vs. NH all)
      ts=ts(kl)*1.0641-0.0764   ; convert to oC wrt 1961-90 (vs. NH land>20)
    end
    2: begin        ; Age-banded MXD
      alltit="Age-banded density NH growing-season reconstruction"
      ; Period to consider
      perst=1402
      peren=1960
;      fac=0.0       ; do not smooth it any further!!!
;      restore,filename='../treeharry/densadj_all(330).idlsave'
;      timey=x
;          ; CONVERSION FACTORS FOR AGE-BANDED MXD, BY REGRESSION ON INSTR.
;      ts=densadj*0.156525   ; converts it from density to temperature anom
      timey=newagetime
      ts=newagets
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
;      ts=ts(kl)-0.140369   ; to convert it oC wrt 1961-90
    end
    3: begin         ; Torn+Yama+Taim
      perst=1
      peren=1993
      openr,2,'tornyamataim.ave'
      readf,2,nnn
      rawdat=fltarr(2,nnn)
      readf,2,rawdat
      close,2
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)*0.1342-0.2761    ; to convert it oC wrt 1961-90
      ;
    end
    4: begin         ; Overpeck
      perst=1600
      peren=1990
      fac=0.18       ; multiply filter by this because record is already 5-yrs
      openr,2,'overpeck_arctic.dat'
      headst=strarr(2)
      readf,2,headst
      readf,2,nnn
      rawdat=fltarr(2,nnn)
      readf,2,rawdat
      close,2
      timey=reverse(reform(rawdat(0,*)))
      ts=reverse(reform(rawdat(1,*)))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)*0.3352-0.0802    ; to convert it oC wrt 1961-90
      ;
    end
    5: begin         ; Crowley & Lowery
      fac=0.8        ; multiply filter by this 'cos record is already smooth
      restore,'/cru/u2/f055/data/paleo/crowley/crowley_lowery.idlsave'
      ; Gets crownyr,crowtimey,crowts
      timey=crowtimey
      ts=crowts
    end
    6: begin         ; Instrumental NH land > 20N, Apr-Sep
      perst=1871
      peren=1997
      openr,2,'../treeharry/nhland20_amjjas.dat'
      readf,2,nnn
      rawdat=fltarr(2,nnn)
      readf,2,rawdat
      close,2
      print,rawdat(0,0)
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
      ;
    end
  endcase
  ;
  ; Now plot the results
  ;
  if fac gt 0. then filter_cru,thalf*fac,/nan,tsin=ts,tslow=tslow $
               else tslow=ts
  ;
  if jdo eq 0 then begin
    plot,timey,tslow,/nodata,$
      xstyle=1,xrange=[1000,2000],xtitle='Year  (AD)',$
      /ystyle,ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
      yrange=[-0.75,0.25]
;      yrange=[-0.61,0.25]
    ttt=newagetime
    tts=newagelow
    tse=newagese
    tkk=where((ttt ge 1402) and (ttt le 1960))
    ttt=ttt(tkk) & tts=tts(tkk) & tse=tse(tkk)
    if doerr eq 1 then begin
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse*2.,reverse(tts+tse*2.)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
      polyfill,xfill,yfill,color=19,noclip=0
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse,reverse(tts+tse)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
      polyfill,xfill,yfill,color=21,noclip=0
    endif
    oplot,!x.crange,[0.,0.],linestyle=1
  endif
  ;
  if ido eq 2 then tslow=newagelow
  oplot,timey,tslow,thick=lthi(ido),color=lcol(ido),linestyle=lsty(ido)
print,min(tslow,/nan)
  ;
  ; For Overpeck, let's interpolate to get annual values from decadal ones
  ;
  if ido eq 4 then begin
    inyr=timey(nyr-1)-timey(0)+1
    itimey=findgen(inyr)+timey(0)
    its=interpol(ts,timey,itimey)
    ts=its
    timey=itimey
  endif
  ;
  ; Store timeseries for later cross-correlation
  ;
  ist=where(alltime eq timey(0))
  alldo(ido,ist(0):ist(0)+n_elements(ts)-1)=ts(*)
  ;
endfor
;
case doerr of
  0: terr=''
  1: terr=' (& errors)'
endcase
legpos=convert_coord([1025],[-0.42],/data,/to_normal)
legord=[0,1,2,5,6]
legtxt=['Jones et al. (1998)','Mann et al. (1999)','Briffa et al. (2000)',$
    'Crowley & Lowery (2000)',$
    'Observations']
legend,position=legpos,legtxt,$
  thick=lthi(legord),color=lcol(legord),linestyle=lsty(legord)
;
end