;
; Synthesises together multiple quasi-NH recons
;
if n_elements(doseas) eq 0 then doseas=1  ; 0=annual, 1=apr-sep, 2=oct-mar
seasname=['annual','aprsep','octmar']
mtit=['ANNUAL (Jan-Dec)','WARM SEASON (Apr-Sep)','COLD SEASON (Oct-Mar)']
dobore=0  ; 0=no boreholes, 1=overlay borehole-based NH estimate
doerr=0   ; 0=no errors, 1=Briffa age errors
dost=1850 ; start year
docorr=0        ; 1=compute correlations  0=don't
thalf=0.       ; filter for high or low-pass
dobigy=1        ; 0=compact yrange  1=big yrange
if doseas ne 1 then dobigy=1.4
doinstrwarm=0   ; 1=use warm-season instrumental series as predictor too!
;
; 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!'
;
lcol=[10,11,12,13,14,16,15,17,50]
if thalf gt 0 then begin
  lthi=[6,6,6,6,6,6,4,6]
endif else begin
  lthi=[2,2,2,2,2,2,4,2]
endelse
;
loadct,39
multi_plot,nrow=2,layout='large'
if !d.name eq 'X' then begin
  lthi=lthi*0.5
  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='brightblue'
def_1color,16,color='orange'
def_1color,17,color='magenta'
def_1color,50,color='mgrey'
;
def_1color,18,color='white'
def_1color,22,color='vlblue'
def_smearcolor,fromto=[18,22]
def_1color,21,color='mdgrey'
def_1color,22,color='lgrey'
;
ndo=8
allnyr=2002
alltime=findgen(allnyr)
alldo=fltarr(ndo,allnyr)*!values.f_nan
lowdo=fltarr(ndo,allnyr)*!values.f_nan
namedo=['Jones ','Mann  ','Briffa','NChron','Overpeck','Crowley & Lowery',$
  'Instr','Esper']
iord=[3,0,4,5,1,7,2,6]
;
for jdo = 0 , ndo-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)
      case doseas of
        0: ts=ts*0.3690-0.2051       ; via regression vs. NH land>20 ann
        1: ts=ts*0.3369-0.1209       ; via regression vs. NH land>20 warm
        2: ts=ts*0.3949-0.2989       ; via regression vs. NH land>20 cold
      endcase
    end
    1: begin        ; Mike Mann's reconstruction (now use his 1000 yr one)
                    ; ALSO, NOW HAVE A LAND-NORTH-OF-20N ANNUAL RECON
      alltit="Mann's Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000
      peren=1980
      doln20=1
      if doln20 eq 0 then begin
        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]
;        ts=ts-0.12            ; convert to oC wrt 1961-90 (vs. NH all)
;        ts=ts*1.0641-0.0764   ; convert to oC wrt 1961-90 (vs. NH land>20)
        case doseas of
          0: ts=ts*1.1106-0.1645       ; via regression vs. NH land>20 ann
          1: ts=ts*0.9532-0.0886       ; via regression vs. NH land>20 warm
          2: ts=ts*1.2834-0.2480       ; via regression vs. NH land>20 cold
        endcase
      endif else begin
        openr,1,'../tree5/mannarea_all.dat'
        nyr=981
        headdat=' '
        rawdat=fltarr(11,nyr)
        readf,1,headdat
        readf,1,rawdat           ;,format='(I6,F11.7)'
        close,1
        timey=reform(rawdat(0,*))
        ts=reform(rawdat(10,*))
        kl=where((timey ge perst) and (timey le peren),nyr)
        timey=timey(kl)
        ts=ts[kl]
;OLD VERSION ts=ts*0.7246-0.0945   ; convert to oC wrt 1961-90 (vs. NH land>20)
        case doseas of
          0: ts=ts*0.8270-0.1771       ; via regression vs. NH land>20 ann
          1: ts=ts*0.6354-0.1061       ; via regression vs. NH land>20 warm
          2: ts=ts*0.9820-0.2602       ; via regression vs. NH land>20 cold
        endcase
      endelse
    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
      case doseas of
        0: ts=ts*0.8996-0.1079       ; via regression vs. NH land>20 ann
        1: ts=ts*0.9281-0.0151       ; via regression vs. NH land>20 warm
        ; this isn't 1.0-0.0 because new instr data modified it
        2: ts=ts*0.7619-0.2271       ; via regression vs. NH land>20 cold
      endcase
    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]
;      ts=ts*0.1342-0.2761    ; to convert it oC wrt 1961-90
      case doseas of
        0: ts=ts*0.1308-0.3650       ; via regression vs. NH land>20 ann
        1: ts=ts*0.1188-0.2664       ; via regression vs. NH land>20 warm
        2: ts=ts*0.1452-0.4746       ; via regression vs. NH land>20 cold
      endcase
    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]
;      ts=ts(kl)*0.3352-0.0802    ; to convert it oC wrt 1961-90
      case doseas of
        0: ts=ts*0.3900-0.1593       ; via regression vs. NH land>20 ann
        1: ts=ts*0.2737-0.0947       ; via regression vs. NH land>20 warm
        2: ts=ts*0.4812-0.2303       ; via regression vs. NH land>20 cold
      endcase
    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.'+$
        seasname[doseas]+'.idlsave'
      ; Gets crownyr,crowtimey,crowts
      timey=crowtimey
      ts=crowts
    end
    6: begin         ; Instrumental NH land > 20N, Apr-Sep (or other season!)
      perst=1881
      peren=1993
      print,'Reading temperatures'
      ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2_18512001.mon.nc')
      tsmon=crunc_rddata(ncid,tst=[1851,0],tend=[2001,11],grid=gtemp)
      ncdf_close,ncid
      nmon=12
      ntemp=gtemp.nt
      nyrtemp=ntemp/nmon
      timey=reform(gtemp.year,nmon,nyrtemp)
      timey=reform(timey(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: ts=mkseason(nhmon,0,11,datathresh=6)
        1: ts=mkseason(nhmon,3,8,datathresh=3)
        2: ts=mkseason(nhmon,9,2,datathresh=3)
      endcase
      warmyr=timey
      warmts=mkseason(nhmon,3,8,datathresh=3)
      case doseas of
        0: warmts=warmts*0.8929-0.1058   ; via regression vs. NH land>20 ann
        1: warmts=warmts*1.0000-0.0000   ; via regression vs. NH land>20 warm
        2: warmts=warmts*0.7365-0.2285   ; via regression vs. NH land>20 cold
      endcase
;      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)
      warmyr=warmyr[kl]
      warmts=warmts[kl]
    end
    7: begin        ; Esper's reconstruction
      alltit="Esper et al.'s Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=831
      peren=1992
      ; Read raw reconstruction
      openr,1,'/cru/u2/f055/data/paleo/esper2002/esper.txt'
      readf,1,nyr
      headst=''
      readf,1,headst
      rawdat=fltarr(7,nyr)
      readf,1,rawdat
      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)
      ; Convert from raw values to deg C relative to 1961-1990
      case doseas of
        0: ts=ts*1.6664-2.1539       ; via regression vs. NH land>20 ann
        1: ts=ts*1.4724-1.8442       ; via regression vs. NH land>20 warm
        2: ts=ts*2.0851-2.7289       ; via regression vs. NH land>20 cold
      endcase
    end
  endcase
  ;
  if fac gt 0. then filter_cru,thalf*fac,/nan,tsin=ts,tslow=tslow $
               else tslow=ts
  ;
  ; 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)
    lowits=interpol(tslow,timey,itimey)
    its=interpol(ts,timey,itimey)
    ts=its
    tslow=lowits
    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(*)
  lowdo(ido,ist(0):ist(0)+n_elements(ts)-1)=tslow(*)
  ;
endfor
;
; Having read in all series, we want to synthesise them together.
;
uselist=[0,1,2,3,4,5,7]     ; don't use 6, the instrumental data
nuse=n_elements(uselist)
;
niter=1     ; 12?
for ij = 0 , niter-1 do begin
;
; Now let's re-calibrate each over the period 1000-1850 (or whatever subset
; each has data for) against the mean series.
;
if ij gt 0 then begin
allnew=fltarr(ndo,allnyr)*!values.f_nan
lownew=fltarr(ndo,allnyr)*!values.f_nan
for jdo = 0 , nuse-1 do begin
  ido=uselist[jdo]
  ts1=reform(lowdo[ido,*])
  ts2=reform(alldo[ido,*])
  kl=where((alltime ge 1600) and (alltime le 1850) and (finite(ts1) ne 0),$
    nkeep)
  print,'re-calibrating'
  acoeff=linfit(ts1[kl],alllow[kl])
  print,namedo[ido],nkeep,acoeff
  lownew[ido,*]=acoeff[0]+acoeff[1]*ts1[*]
  allnew[ido,*]=acoeff[0]+acoeff[1]*ts2[*]
endfor
lowdo[uselist,*]=lownew[uselist,*]
alldo[uselist,*]=allnew[uselist,*]
endif
;
; Compute their average and calibrate this against the instrumental data,
; estimating uncertainty/errors too!
;
allts=alldo[uselist,*]
alltarg=reform(alldo[6,*])
allave=multits_calib(allts,alltime,alltarg,alltime,calibper=[1881,1960],$
  stderr=serr,no_anomaly=0,no_varadj=0,smoothlist=[thalf],plotdiag=0)
;
; Optionally keep the uncalibrated average of the multiple reconstructions
;
if ij eq 0 then begin
  print,'keeping average'
  keepave=total(allts,1,/nan)/float(total(finite(allts),1))
  filter_cru,thalf,/nan,tsin=keepave,tslow=keeplow
endif
;
; Make smoothed series for plotting
;
filter_cru,thalf,/nan,tsin=allave,tslow=alllow
!p.multi[0]=0
;
;print,'Calibrate average vs. observed:',r,acoeff
  ;
  ; Now plot the results
  ;
if (ij eq 0) or (ij eq niter-1) then begin
for jdo = 0 , ndo-1 do begin
  ido=iord(jdo)
  tslow=reform(lowdo[ido,*])
  ;
  if jdo eq 0 then begin
    yrange=[-0.75-0.2*dobore-(doseas eq 2)*0.2,0.4+dobigy*0.05+$
        (doseas eq 2)*0.1]
    if thalf eq 100. then yrange=[-0.95,0.2]
    yrange=[-0.8,0.6]
;    yrange=[-0.75,0.25]
    pause
    plot,alltime,tslow,/nodata,$
      xstyle=1,xrange=[dost,2001],xtitle='Year  (AD)',$
      /ystyle,ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
      yrange=yrange
;    ttt=newagetime
;    tts=newagelow
;    tse=newagese
;    tkk=where((ttt ge 1402) and (ttt le 1960))
;    ttt=ttt(tkk) & tts=tts(tkk) & tse=tse(tkk)
    ttt=alltime
    tts=alllow
    tse=serr
    if doerr eq 1 then begin
print,'PLOTTING ERRORS!'
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse*2.,reverse(tts+tse*2.)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
kerr=where(xfill ge 1600)
      polyfill,xfill[kerr],yfill[kerr],color=21,noclip=0
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse,reverse(tts+tse)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
kerr=where(xfill ge 1600)
      polyfill,xfill[kerr],yfill[kerr],color=22,noclip=0
    endif
    oplot,!x.crange,[0.,0.],linestyle=1
;
; Optionally overlay pre-computed borehole temperature estimates
;
if dobore eq 1 then begin
  restore,'/cru/u2/f055/data/paleo/borehole/borehole_NH.idlsave'
  ; Gets: boretime, boresite, boregrid
  ; Compute and apply an offset so that the borehole line crosses zero
  ; in 1975.5 (i.e., mid-way through the 1961-1990 reference period)
  v=boresite[0:1]
  voff=(v[0]-v[1])*(2000.-1975.5)/(2000.-1900.)-v[0]
;  oplot,boretime,boresite+voff,color=50,thick=6,linestyle=2
  v=boregrid[0:1]
  voff=(v[0]-v[1])*(2000.-1975.5)/(2000.-1900.)-v[0]
  oplot,boretime,boregrid+voff,color=50,thick=6,noclip=1
endif
  ;
  endif
  ;
  oplot,alltime,tslow,thick=lthi(ido),color=lcol(ido)
endfor
;
if doinstrwarm eq 1 then begin
  filter_cru,thalf,/nan,tsin=warmts,tslow=warmlow
  oplot,warmyr,warmlow,thick=4
endif
;
case doerr of
  0: terr=''
  1: terr=' (& errors)'
endcase
case dost of
  1: iiipos=50
  850: iiipos=880
  950: iiipos=950
  1000: iiipos=1025
  1400: iiipos=1415
  1800: iiipos=1940
  else: iiipos=dost
endcase
if dost le 1000 then jjjpos=-0.39 else jjjpos=0.35
iiipos=!x.crange[0]
jjjpos=!y.crange[1]
legpos=convert_coord([iiipos],[jjjpos],/data,/to_normal)
legord=[4,0,1,5,3,2,7,6]
;legtxt=['Jones et al. (1998)','Mann et al. (1999)','This study'+terr,$
;    'Briffa (2000)','Overpeck et al. (1997)','Crowley & Lowery (2000)',$
;    'Observations']
;legtxt=['Jones et al. (1998)','Mann et al. (1999)','Briffa et al. (2000)',$
;    'Briffa (2000)','Overpeck et al. (1997)','Crowley & Lowery (2000)',$
;    'Observations','Esper et al. (2002)']
;legend,position=legpos,legtxt(legord),$
;  thick=lthi(legord),color=lcol(legord)
legtxt=['Jones98','Mann99','Briffa01',$
    'Briffa00','Overpeck97','Crowley00',$
    'Obs','Esper02']
legend,position=legpos,legtxt(legord),horiz=1,$
  margin=0.5,pspacing=1,thick=lthi(legord),color=lcol(legord)
xyouts,dost+75,!y.crange[1]-0.17,mtit[doseas]
;
if (ij ne niter-1) or (ij eq 0) then begin
;  oplot,alltime,alllow,thick=4+4*(thalf gt 0)
;  oplot,alltime,keeplow,thick=2+(thalf gt 0)
endif
endif
;
endfor
;
end