;
; Compares the calibrated regional series of Hugershoff and Age-Banded
; MXD records
;
showband=0       ; 0=do Hug & Hug errors, 1=do Hug, Age & Age errors, -1=Age
do9=3             ; 0=do 7 regions & ALL,  1=do 9 regions, 2=do 5 regions
dolongeuro=0      ; 0=don't plot lon euro instrumental series, 1=do
cut1960=0         ; 0=plot full data, 1=cut off at 1960
docol=1           ; 0=B&W 1=color
;
if (showband ne 0) and (do9 eq 0) then message,'Age-banded ALL is not PCR vers'
;
if showband eq 0 then begin
  maints='reg'
  plusts='band'
endif else begin
  maints='band'
  plusts='reg'
endelse
;
thalfwant=10.
addst='sm'+string(thalfwant,format='(I2.2)')
;
; Read in long instrumental records from Europe
;
restore,$
filename='/cru/u2/f055/data/obs/station/surface/longeuro/longeuro_seas.idlsave'
;  Gets: nyr,timey,nseas,seasname,nstat,statname,nreg,regname,seasts,regseas
;
liyear=timey
lineur=reform(regseas(0,*,5))
liseur=reform(regseas(2,*,5))
filter_cru,thalfwant,/nan,tsin=lineur,tslow=neurlow
filter_cru,thalfwant,/nan,tsin=liseur,tslow=seurlow
ml=where(liyear lt 1750.)
;neurlow(ml)=!values.f_nan
;seurlow(ml)=!values.f_nan
;
; Read in the age-banded calibrated series ERRORS
;
;addst=''
;restore,filename=maints+'temp'+addst+'_mxd_errors.idlsave'
;; Gets: dosmooth,thalf,calregse,regname,timey
;banderryr=timey
;
; Check that the errors were run for the same filter length
;
;if dosmooth eq 0 then message,'Need smoothed errors!'
;if thalf ne thalfwant then message,'Need errors for same smoothing!'
;
; Read in the age-banded calibrated series
;
restore,filename=maints+'temp_calibrated.allversion.idlsave'
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
calbandts=calregts
bandyr=timey
;if n_elements(banderryr) ne n_elements(bandyr) then message,'Oh dear'
;
; Read in the Hugershoff calibrated series
;
restore,filename=plusts+'temp_calibrated.idlsave'
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
hugyr=timey
;
; Prepare for plotting
;
loadct,39
if do9 eq 3 then nrrr=4 else nrrr=9
multi_plot,nrow=nrrr,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=20
endelse
!y.margin=[0,0]
!y.omargin=[4,4]
!x.ticklen=0.08
if docol eq 1 then begin
  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'
endif else begin
  def_1color,20,color='mgrey'
  def_1color,21,color='palegrey'
  def_1color,22,color='black'
  def_1color,23,color='black'
  def_1color,24,color='black'
endelse
;
; Repeat for each region
;
cutlist=[1000,1620,1000,1460,1650,1000,1580,1710,1740,1000]
partno='('+['a','b','c','d','e','f','g','h','i','j']+') '
case do9 of
  0: dolist=[0,1,2,3,6,7,8,9]
  1: dolist=indgen(9)
  2: dolist=[0,2,3,7,8]
  3: dolist=[9]
endcase
ndo=n_elements(dolist)
for ido = 0 , ndo-1 do begin
  ireg=dolist(ido)
  ;
  if do9 eq 2 then begin
    if ido ne 2 then ytit='' $
                else ytit='Temperature anomaly wrt. 1961-90 (!Uo!NC)'
  endif else begin
    if ((ido mod 2) eq 0) then ytit='' $
                          else ytit='Temperature anomaly wrt. 1961-90 (!Uo!NC)'
  endelse
                ytit='Temperature anomaly wrt. 1961-90 (!Uo!NC)'
  ;
  if ireg eq nreg-1 then begin
    yr=[-1.7,0.7]
    ypos=-1.6
;    !p.multi=[1,1,5,0,0]
  endif else begin
    yr=[-2.2,1.0]
    ypos=-1.8
  endelse
  ;
  if cut1960 eq 0 then maxx=1994 else maxx=1960
  if do9 eq 3 then begin
    minx=1650
    ntic=7
    yr=[-2.5,1.5]
  endif else begin
    minx=1400
    ntic=12
    yr=[-1.99,0.99]
  endelse
  minx=1870
  maxx=1995
  yr=[-1.,0.9]
  plot,[0,1],/nodata,$
    /xstyle,xrange=[minx,maxx],xtickformat='nolabels',$
;    xticks=ntic-1,xtickv=findgen(ntic)*50+minx,xminor=5,$
    /ystyle,yrange=yr,$
    ytitle=ytit
  ;
  serr=calregse(*,ireg)
  ;
  tem1=reform(tempregts(*,ireg))
  filter_cru,thalf,/nan,tsin=tem1,tslow=temlow
  ;
  age1=reform(calbandts(*,ireg))
  filter_cru,thalf,/nan,tsin=age1,tslow=agelow
  ;
  hug1=reform(calregts(*,ireg))
  filter_cru,thalf,/nan,tsin=hug1,tslow=huglow
  ;
  ml=where((finite(age1) eq 0) and (bandyr le 1800),nmiss)
  if nmiss gt 0 then begin
    lastmiss=ml(nmiss-1)
    print,'Removing data up to',bandyr(lastmiss)
    age1(0:lastmiss)=!values.f_nan
  endif
  ml=where(finite(age1) eq 0,nmiss)
  if nmiss gt 0 then begin
    agelow(ml)=!values.f_nan
    huglow(ml)=!values.f_nan
  endif
  ;
;  xfill=[bandyr,reverse(bandyr)]
;  yfill=[agelow-serr*2.,reverse(agelow+serr*2.)]
;  kl=where(finite(yfill),nkeep)
;  yfill=yfill(kl)
;  xfill=xfill(kl)
;  polyfill,xfill,yfill,color=20,noclip=0
  ;
;  xfill=[bandyr,reverse(bandyr)]
;  yfill=[agelow-serr,reverse(agelow+serr)]
;  kl=where(finite(yfill),nkeep)
;  yfill=yfill(kl)
;  xfill=xfill(kl)
;  polyfill,xfill,yfill,color=21,noclip=0
  ;
  axis,yaxis=0,/ystyle,ytickformat='nolabels'
  axis,yaxis=1,/ystyle,ytickformat='nolabels'
  ;
;  if !d.name eq 'PS' then device,font_size=9
; xyouts,1910,ypos,partno(ido)+regname(ireg)
;  if !d.name eq 'PS' then device,font_size=12
  ;
  if ido eq 0 then axis,xaxis=1,xtitle='Year',/xstyle
  if ido eq ndo-1 then axis,xaxis=0,xtitle='Year',/xstyle
  ;
  oplot,!x.crange,[0.,0.],linestyle=1
  ;
  oplot,temptimey,tem1,thick=3+(1-docol),color=22
  if showband gt 0 then oplot,hugyr,huglow,thick=1+2*docol,color=23
  oplot,bandyr,age1,thick=2+docol
kkk1=where((temptimey ge 1881) and (temptimey le 1990))
kkk2=where((bandyr ge 1881) and (bandyr le 1990))
print,correlate(tem1[kkk1],age1[kkk2])
  if dolongeuro ne 0 then begin
    if ireg eq 0 then oplot,liyear,neurlow,thick=3,color=24,linestyle=2-2*docol
    if ireg eq 1 then oplot,liyear,seurlow,thick=3,color=24,linestyle=2-2*docol
  endif
  ;
  if showband ne 0 then oplot,replicate(cutlist(ireg),2),!y.crange,thick=3,color=22
  ;
endfor
;
;!p.multi=[1,1,6,0,0]
;
;plot,[0,1],/nodata,$
;  /xstyle,xrange=[1400,1994],xtickformat='nolabels',$
;  /ystyle,yrange=[-2.2,1.0],$
;  ytitle=ytit
;
end