;
; Compares the calibrated regional series of Age-Banded
; MXD records, for inverse and classical regression.
;
ploton,1
;
showband=1       ; 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
                  ; 3=do 1 region & bigger, 4=do 2 regions & bigger
                  ; 5=do 7 regions
dolongeuro=1      ; 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)+regseas(1,*,5))/2.
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
;
if maints eq 'band' then addst=''
fn1=maints+'temp'+addst+'_mxd_errors.idlsave'
print,fn1
restore,filename=fn1
; 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
;
;if maints eq 'band' then begin
;  addver=''
;endif else begin
;  if regname[9] eq 'ALL' then addver='allversion.' else addver='nhversion.'
;endelse
;restore,filename=maints+'temp_calibrated.'+addver+'idlsave'
restore,filename='bandtemp_calibrated.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
;
;if plusts eq 'band' then begin
;  addver=''
;endif else begin
;  if regname[9] eq 'ALL' then addver='allversion.' else addver='nhversion.'
;endelse
;restore,filename=plusts+'temp_calibrated.'+addver+'idlsave'
restore,filename='bandtemp_classiccalibrated.idlsave'
; Gets: nyr,nreg,calregts,regname,timey,tempregts,tempnyr,temptimey
hugyr=timey
;
; Prepare for plotting
;
loadct,39
if do9 ge 3 then nrrr=3 else nrrr=9
if do9 eq 5 then nrrr=7
if do9 eq 4 then fs=18 else fs=12
if do9 eq 3 then fs=20
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=fs
endelse
!y.margin=[0,0]
!y.omargin=[4,4]
!x.ticklen=0.08
if docol eq 1 then begin
  def_1color,20,color='mlorange'
  def_1color,21,color='yellow'
  def_1color,22,color='magenta'
  def_1color,23,color='brightblue'
  def_1color,24,color='dgreen'
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]
  4: dolist=[0,1]
  5: dolist=[0,1,2,3,6,7,8]
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
  ;
  if ireg eq nreg-1 then begin
    yr=[-1.7,0.7]
    ypos=-1.6
    if do9 ne 3 then !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
    minx=1800
    ntic=7
    yr=[-2.2,0.8]
  endif else begin
    if do9 eq 4 then begin
      minx=1500
      ntic=6
      yr=[-2.2,1.2]
    endif else begin
      minx=1400
      ntic=12
      yr=[-2.99,1.49]
    endelse
  endelse
;minx=1700
;yr=[-1.7,1.2]
print,maxx
  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
  ;
  kcut1=where(hugyr ge cutlist[ireg])
  kcut2=where(bandyr ge cutlist[ireg])
  ;
  xfill=[bandyr[kcut2],reverse(bandyr[kcut2])]
  yfill=[agelow[kcut2]-serr[kcut2]*2.,reverse(agelow[kcut2]+serr[kcut2]*2.)]
  kl=where(finite(yfill),nkeep)
  yfill=yfill(kl)
  xfill=xfill(kl)
  polyfill,xfill,yfill,color=20,noclip=0
  ;
  xfill=[bandyr[kcut2],reverse(bandyr[kcut2])]
  yfill=[agelow[kcut2]-serr[kcut2],reverse(agelow[kcut2]+serr[kcut2])]
  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=fs-3
  xyouts,1910,ypos,partno(ido)+regname(ireg)
  if !d.name eq 'PS' then device,font_size=fs
  ;
  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
  ;
  th=4+(1-docol)
  oplot,temptimey,temlow,thick=th+1
  oplot,temptimey,temlow,thick=th,color=22
  if showband gt 0 then oplot,hugyr[kcut1],huglow[kcut1],thick=2+2*docol,color=23,noclip=1
  oplot,bandyr[kcut2],agelow[kcut2],thick=5+docol
  if dolongeuro ne 0 then begin
    if ireg eq 0 then oplot,liyear,neurlow,thick=6,color=24,linestyle=2-2*docol
    if ireg eq 1 then oplot,liyear,seurlow,thick=6,color=24,linestyle=2-2*docol
  endif
  ;
  ;if showband ne 0 then oplot,replicate(cutlist(ireg),2),!y.crange,thick=3,color=22
  ;
  print,regname(ireg)
  filter_cru,25.,/nan,tsin=age1,tshigh=a25
  filter_cru,25.,/nan,tsin=hug1,tshigh=h25
  nele=min([n_elements(h25),n_elements(a25)])
  h25=h25(0:nele-1)
  a25=a25(0:nele-1)
  bandyr=bandyr[0:nele-1]
  rrr=mkcorrelation(a25,h25,bandyr)
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1881,1960])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1400,1880])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1400,1499])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1500,1599])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1600,1699])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1700,1799])]
  rrr=[rrr,mkcorrelation(a25,h25,bandyr,refperiod=[1800,1899])]
  print,'Hi-pass r(age,hug)=',string(rrr,format='(8F8.2)')
  ;
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
;
plotoff
spawn,'ps2web'
;
end