;
; Computes lists of grid-boxes in each region
; BUT MAKES IT TIME-VARYING, ACCORDING TO WHEN THE DENSITY CHRONOLOGIES
; ARE MISSING/NON-MISSING.
;
;----------------------------------------------------
;
; Get chronology locations and timeseries (to tell when they are missing!)
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_varget,ncid,'longitude',xlon
ncdf_varget,ncid,'latitude',ylat
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_close,ncid
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
; Get region info
;
restore,'reglists.idlsave'
;
; Get IPCC instrumental dataset grid and define arrays
;
g=def_grid(/ipcc)
maxbox=99           ; alter this if changing regions/chronologies
maxtime=1992-1881+1 ; to keep space to a minimum, only do the instrumental yrs
boxlists=fltarr(maxbox,maxtime,nreg)*!values.f_nan
boxcount=fltarr(maxtime,nreg)
;
; For each region, locate the grid-boxes with data in
;
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)
  allt=treelist(0:ntree(ireg)-1,ireg)     ; subset of tree indices for region
  for i = 0 , maxtime-1 do begin
    maskyr=density(i+1881-1600,allt)      ; corresponding densities for year i
    keeplist=where(finite(maskyr),nkeep)  ; which of the subset are there?
    print,i+1881,nkeep,format='($,2I5)'
    if nkeep eq 0 then begin
      boxcount(i,ireg)=0.
    endif else begin
      gott=allt(keeplist)                 ; the available subset of the subset!
      statx=xlon(gott)
      staty=ylat(gott)
      statval=fltarr(nkeep)+1.
      fdout=gridit(g.nx,g.ny,g.x,g.y,statx,staty,statval)
      boxele=where(finite(fdout),ngot)
      boxlists(0:ngot-1,i,ireg)=boxele(*)
      boxcount(i,ireg)=ngot
    endelse
  endfor
  print
endfor
;
save,filename='timevary.idlsave',nreg,boxlists,maxbox,maxtime,boxcount,$
  regname
;
end