;
; Computes rbar for each region, for instrumental
;
; Define dimensions
;
yrst=1856
yren=1992
nyr=yren-yrst+1
nmon=12
g=def_grid(/ipcc)
;
; Get masks of boxes to use for each region
;
restore,filename='regboxes.idlsave'
;
; Convert to a list of grid-points
;
ntree=intarr(nreg)
for i = 0 , nreg-1 do begin
  fdmask=reform(boxlists(*,*,i),g.nx*g.ny)
  dummy=where(finite(fdmask),ndummy)
  ntree(i)=ndummy
endfor
maxbox=max(ntree)
densmon=fltarr(nyr,nmon,maxbox,nreg)
densann=fltarr(nyr,maxbox,nreg)
treelist=fltarr(maxbox,nreg)
for i = 0 , nreg-1 do begin
  fdmask=reform(boxlists(*,*,i),g.nx*g.ny)
  dummy=where(finite(fdmask),ndummy)
  treelist(0:ntree(i)-1,i)=dummy(*)
endfor
;
; Read in monthly data and keep the groups of grid boxes only
;
openr,1,'/cru/dave3/f055/detect/obsdat/tair_monthly.dat'
for iyr = 0 , nyr-1 do begin
  for imon = 0 , nmon-1 do begin
    fd=rdtair(72,36,imon=im,iyr=iy,lun=1)
    fd=reform(fd,g.nx*g.ny)
    print,iyr,iy,imon,im
    for i = 0 , nreg-1 do begin
      densmon(iyr,imon,0:ntree(i)-1,i)=fd(treelist(0:ntree(i)-1,i))
    endfor
  endfor
endfor
close,1
;
; Convert to a long-season mean (just let it create divide by zero errors)
; BUT DO MONTHLY RBAR AS WELL!
;
densamjjas = total(densmon(*,3:8,*,*),2,/nan) / $
  float(total(finite(densmon(*,3:8,*,*)),2,/nan))
densmjja = total(densmon(*,4:7,*,*),2,/nan) / $
  float(total(finite(densmon(*,4:7,*,*)),2,/nan))
densjj = total(densmon(*,5:6,*,*),2,/nan) / $
  float(total(finite(densmon(*,5:6,*,*)),2,/nan))
;
; Repeat for each region
;
amjjasrbar=fltarr(nreg)
mjjarbar=fltarr(nreg)
jjrbar=fltarr(nreg)
monrbar=fltarr(nreg,nmon)
for i = 0 , nreg-1 do begin
  ;
  ; Extract number of chronologies and the timeseries
  ;
  print,'Doing region ',regname(i),' with ntree=',ntree(i)
  nsize=ntree(i)
  ;
  inx=reform(densamjjas(*,0:nsize-1,i))
  amjjasrbar(i)=rbar(transpose(inx))
  ;
  inx=reform(densmjja(*,0:nsize-1,i))
  mjjarbar(i)=rbar(transpose(inx))
  ;
  inx=reform(densjj(*,0:nsize-1,i))
  jjrbar(i)=rbar(transpose(inx))
  ;
  for imon = 0 , nmon-1 do begin
    inx = reform(densmon(*,imon,0:nsize-1,i))
    monrbar(i,imon)=rbar(transpose(inx))
  endfor
  ;
endfor
;
save,filename='rbar_monthly.idlsave',nreg,regname,amjjasrbar,jjrbar,$
  monrbar,mjjarbar
;
end