;
; Computes rbar for each region (not composited regions) for densities, and
; also saves the grand correlation matrices
;
; Get MXD 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 to use only a common period
;
keeplist=where((x ge 1881) and (x le 1960),nkeep)
if nkeep gt 0 then begin
  density=density(keeplist,*)
  x=x(keeplist)
endif
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
; Get lists of trees in each region
;
restore,filename='reglists.idlsave'
;
; Repeat for each region
;
for i = 3 , nreg-1 do begin      ; ignore 0, 1, and 2 (ALL, NORTH & SOUTH)
  ;
  ; 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))
  ;
  ; Define a correlation matrix
  ;
  grandr=fltarr(nsize,nsize)
  ;
  ; Do each tree against all others (including itself)
  ;
  for j = 0 , nsize-1 do begin
    print,j,format='($,I5)'
    iny=reform(inx(*,j))
    rowr=manycorr(inx,iny,/nan,nval=nval)     ; compute r values
    grandr(*,j)=rowr(*)
  endfor
  print
  ;
  save,filename='rbartd_dens_'+regname(i)+'.idlsave',grandr
  ;
endfor
;
end