;
; 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=['allsites']
nreg=n_elements(regname)
;
; Specify the number of bands used in each region
;
nband=intarr(nreg)
ncol=intarr(nreg)
;
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,ncol=2,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=12
endelse
!p.multi(4)=1
!y.margin=[0,0]
!y.omargin=[4,2]
!x.margin=[0,0]
!x.omargin=[10,4]
;
; Repeat for each region
;
for ireg = 0 , nreg-1 do begin
  ;
  ; Read raw data
  ;
  fn='mxd.'+regname(ireg)+'.BS5T1S1.17002000.sa.dat1950'
  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
  ;
  ; 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.
  ;
  pause
  plot,timey,timey,/nodata,$
    /xstyle,/ystyle,yrange=[-41,3],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
  ;
  ; 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