;
; Computes rbar for ALL region only, for density
;
; Get ring width timeseries and locations
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'fraction',weight
ncdf_close,ncid
;
; Reduce length of it
;
help,density,x
keeplist=where((x ge 1881) and (x le 1960),nkeep)
if nkeep gt 0 then begin
  density=density(keeplist,*)
  x=x(keeplist)
endif
help,density,x
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
; Get lists of trees in each region
;
restore,filename='reglists.idlsave'
;
; Define arrays to hold the mean rbar and the grandr
; arrays for each region
;
nreg=1
nmaxr=max(ntree)
grandr=fltarr(nmaxr,nmaxr,nreg)*!values.f_nan
allrbar=fltarr(nreg)
;
; Repeat for each region
;
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=density(*,treelist(0:nsize-1,i))
  ;
  ; Do each tree against all others (including itself)
  ;
  for j = 0 , nsize-1 do begin
    print,j
    iny=reform(inx(*,j))
    rowr=manycorr(inx,iny,/nan,nval=nval)     ; compute r values
    grandr(0:nsize-1,j,i)=rowr(*)
  endfor
  ;
  ; Total it (for nsize by nsize) and average it, but only do i<>j.
  ;
  totr=0.
  totn=0.
  for j = 1 , nsize-1 do begin
    for k = 0 , j-1 do begin
      if finite(grandr(j,k)) then begin
        totr=totr+grandr(j,k)
        totn=totn+1.
        if grandr(j,k) eq 1. then message,'Ooops!'
      endif
    endfor
  endfor
  totr=totr/totn
  ;
  allrbar(i)=totr
  ;
endfor
;
grandr=reform(grandr)
allrbar=reform(allrbar)
save,filename='rbar_dens_18811960.idlsave',nreg,regname,allrbar,grandr
;
end