;
; Plots all bands from all regions, low-pass filtered.
;
; Define region names.  Names are different to those I use, but order
; should be the same to make later comparison easier.
;
regname=['northwestern.europe',$
         'southwestern.europe',$
         'northern.siberia',$
         'eastern.siberia',$
         'centralsiberia',$
         'tibetianplateau',$
         'western.northamerica',$
         'northwestern.canada',$
         'centraleastern.canada']
nreg=n_elements(regname)
nreg=1
;
; Specify the number of bands used in each region
;
; **These values are now overwritten, because they're lower for the early
; **fixed grids!
nband=[19,20,21,24,18,20,22,17,17]
ncol=nband*3 + 1
;
allnyr=2000
allyr=findgen(allnyr)+1.
kper=where((allyr ge 1400) and (allyr le 1995),nyr)
timey=allyr(kper)
;
loadct,39
multi_plot,nrow=1,layout='centred'
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=10
endelse
def_1color,20,color='mdgrey'
;
; Repeat for each region
;
for ireg = 0 , nreg-1 do begin
  ;
  ; Read raw data
  ;
  fn='mxd.'+regname(ireg)+'.all.BS5T1S1.17002000.sa.dat'
  print,fn
  ;
  ; Check the number of columns
  ;
  spawn,'head -1 '+fn+' | wc -w',osoutput
  ncol(ireg)=fix(osoutput)
  nband(ireg)=(ncol(ireg)-1)/3
  ;
  openr,1,fn
  headlin=fltarr(ncol(ireg))
  readf,1,headlin
  rawdat=fltarr(ncol(ireg),allnyr)
  readf,1,rawdat
  close,1
  ;
  headlin=reform(headlin(1:ncol(ireg)-1),3,nband(ireg))
  bndst=string(reform(headlin(0,*)),format='(I3)')
  bnden=string(reform(headlin(1,*)),format='(I3)')
  bndname=bndst+'-'+bnden
  ytv=-findgen(nband(ireg))*2.
  ;
  ; Just keep the 1400 to 1994 portions
  ;
  rawdat=rawdat(*,kper)
  ;
  ; Locate missing data
  ;
  ml=where(rawdat eq 9990.,nmiss)
  rawdat(ml)=!values.f_nan
  ;
  collist=indgen(nband(ireg))*3+1
  tsdat=rawdat(collist,*)      ; values
  ntree=rawdat(collist+2,*)    ; number of trees
  ;
  tottree=total(ntree,1)*0.01-40.
  ;
  plot,timey,timey,/nodata,$
    /xstyle,xtitle='Year',xmargin=[10,8],$
    /ystyle,yrange=[-41,3],yticks=nband(ireg)-1,ytickv=ytv,$
    ytickname=bndname,ytitle='Maximum latewood density averaged in age bands'
;   title=regname(ireg)
  ;
  nnn=nband(ireg) < 20
  for k = 0 , nnn-1 do begin
    ;
    yval=reform(tsdat(k,*))
    ynum=reform(ntree(k,*))
    misslist=where(ynum lt 2,nmiss)
    y2=yval
    if nmiss gt 0 then y2(misslist)=!values.f_nan
    filter_cru,10.,/nan,tsin=yval,tslow=ylow
    filter_cru,10.,/nan,tsin=y2,tslow=y2low
    ;
    fac=float(k)*(-2.)
    oplot,timey,ylow+fac
    oplot,timey,y2low+fac,thick=3
    oplot,!x.crange,fac+[0.,0.],linestyle=1
    ;
  endfor
  oplot,timey,tottree,thick=4,color=20
  axis,yaxis=0,yticks=4,ytickv=-40+indgen(5),ytickformat='nolabels'
  axis,yaxis=1,yticks=4,ytickv=-40+indgen(5),$
    ytickname=string(findgen(5)*100,format='(I3)')
  xyouts,2050,-38,align=0.5,orient=90.,'No. of tree cores'
  oplot,!x.crange,fac+[0.,0.],linestyle=1
  ;
  ; Compute correlations between bands
  ;
;  rbar=mkrbar(tsdat(0:nnn-1,*),grandr=corrmat)
;  fff='('+string(nnn,format='(I2)')+'I4)'
;  print,corrmat*100.,format=fff
  ;
  ; Compute correlations between each band and the 100-300 yr mean timeseries
  ;
;  refts=tsdat(10:13,*)
;  reftot=total(refts,1,/nan)
;  refnum=total(finite(refts),1)
;  refts=reftot/float(refnum)
  ;
;  for k = 0 , nnn-1 do begin
;    print,mkcorrelation(refts,reform(tsdat(k,*)))*100,format='($,I4)'
;  endfor
;  print
  ;
endfor
;
end