;
; Combines the multi-PCR models into one reconstruction, including
; combining and interpolating the time-dependent errors.
; Plots the timescale dependent errors (tsde), and overlays Hugershoff curves.
;
; Prepare for plotting
;
loadct,39
multi_plot,nrow=4,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,bold=0,font_size=14
endelse
def_1color,20,color='lsand'
def_1color,21,color='lyellow'
def_1color,22,color='red'
def_1color,23,color='mblue'
def_1color,24,color='lgreen'
;
; Read in the Hugershoff calibrated series
;
restore,filename='regtemp_calibrated.idlsave'
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
hugyr=timey
;
for jjj = 0 , 1 do begin
  case jjj of
    0: begin
      dosmooth=0
      thalf=0
      xr=[1871,1994]
      yr=[-1.5,1.]
      end
    1: begin
      dosmooth=1
      thalf=10
      xr=[1400,1994]
      yr=[-2.2,1.]
      end
    2: begin
      dosmooth=1
      thalf=50
      end
   endcase
;
regname='ALL'       ; 'ALL' or 'NH'
;dosmooth=1          ; 0=errors for interannual data 1=for smoothed data
;thalf=25.
;
fnhead='bandtemp'+regname
if dosmooth eq 1 then smst='sm'+string(thalf,format='(I2.2)') else smst=''
fntail=smst+'_calpcr.idlsave'
;
; Read in the different PCR models for different fixed grids
;
perlist=[   0,   0,   0,   1,   2,   3,   4,   5,   5,   6,   6,   7,   8,   9]
fixlist=[1950,1800,1700,1700,1700,1700,1600,1600,1500,1500,1400,1950,1988,1988]
nver=n_elements(perlist)
;
for jver = 0 , nver-1 do begin
  iver=perlist(jver)
  verst=string(iver,format='(I1)')
  ifix=fixlist(jver)
  fixst=string(ifix,format='(I4)')
  ;
  fn=fnhead+verst+'_fixed'+fixst+fntail
  print,fn
  restore,filename=fn
  ; Gets: nyr,nhtit,yrmxd,prednh,fullnh,predserr,nhtemp,tempyr
  ;
  if jver eq 0 then begin
    allnyr=nyr
    alltimey=yrmxd
    allts=fltarr(nyr,nver)
    allse=fltarr(nyr,nver)
  endif
  allts(*,jver)=prednh(*)
  allse(*,jver)=predserr(*)
  ;
endfor
;
; Make a combined reconstruction using different models for different periods,
; and different fixed grid error estimates
;
prednh=fltarr(nyr) & prednh(*)=!values.f_nan
predse=fltarr(nyr) & predse(*)=!values.f_nan
;
for iyr = 0 , nyr-1 do begin
  jyr=alltimey(iyr)
  ;
  ; Identify which model to use
  ;
  iuse=13
  if jyr le 1991 then iuse=12
  if jyr le 1988 then iuse=11
  if jyr le 1981 then iuse=0
  if jyr le 1682 then iuse=3
  if jyr le 1666 then iuse=4
  if jyr le 1624 then iuse=5
  if jyr le 1603 then iuse=6
  if jyr le 1587 then iuse=7
  if jyr le 1479 then iuse=9
  ;
  ; Use the selected model prediction for this year
  ;
  if jyr ge 1402 then begin
    prednh(iyr)=allts(iyr,iuse)
    ;
    ; Also use the selected model to provide the standard error, UNLESS we
    ; are in one of the following periods
    ;
    se=allse(iyr,iuse)
    if (jyr ge 1683) and (jyr le 1700) then begin
      se=allse(iyr,2)
    endif
    if (jyr ge 1701) and (jyr le 1800) then begin
      se=interpol([allse(iyr,2),allse(iyr,1)],[1700,1800],[jyr])
    endif
    if (jyr ge 1801) and (jyr le 1899) then begin
      se=interpol([allse(iyr,1),allse(iyr,0)],[1800,1900],[jyr])
    endif
    if (jyr ge 1480) and (jyr le 1499) then begin
      se=allse(iyr,8)
    endif
    if (jyr ge 1500) and (jyr le 1587) then begin
      se=interpol([allse(iyr,8),allse(iyr,7)],[1500,1600],[jyr])
    endif
    if (jyr ge 1402) and (jyr le 1479) then begin
      se=interpol([allse(iyr,10),allse(iyr,9)],[1400,1500],[jyr])
    endif
    predse(iyr)=se
    ;
  endif
  ;
endfor
;
plot,[0,1],/nodata,$
  /xstyle,xrange=xr,xtitle='Year',$
  /ystyle,yrange=yr,ytitle='Temperature (!Uo!NC wrt 1961-90)';,$
;  title=regname+'  '+smst
;
if dosmooth eq 1 then begin
  filter_cru,thalf,/nan,tsin=nhtemp,tslow=temlow
  filter_cru,thalf,/nan,tsin=prednh,tslow=agelow
endif else begin
  agelow=prednh
  temlow=nhtemp
endelse
;
timey=alltimey
serr=predse
xfill=[timey,reverse(timey)]
yfill=[agelow-serr*2.,reverse(agelow+serr*2.)]
kl=where(finite(yfill),nkeep)
yfill=yfill(kl)
xfill=xfill(kl)
if jjj gt 0 then polyfill,xfill,yfill,color=20,noclip=0
;
xfill=[timey,reverse(timey)]
yfill=[agelow-serr,reverse(agelow+serr)]
kl=where(finite(yfill),nkeep)
yfill=yfill(kl)
xfill=xfill(kl)
if jjj gt 0 then polyfill,xfill,yfill,color=21,noclip=0
;
oplot,!x.crange,[0.,0.],linestyle=1
;
oplot,tempyr,temlow,thick=2+3.*dosmooth,color=22
oplot,timey,agelow,thick=2+3.*dosmooth
;
;if dosmooth eq 1 then begin
  hugts=reform(calregts(*,9))
  if dosmooth eq 1 then filter_cru,thalf,/nan,tsin=hugts,tslow=huglow $
                   else huglow=hugts
  oplot,hugyr,huglow,thick=2+3.*dosmooth,color=23
;endif
;
endfor
;
end