;
; Computes lists of grid-boxes in each region, and
; lists of grid-boxes in each region that also have
; MXD sites in them.
;
;***********
; boxregs=boxes in each big region
; boxprec=boxes in each big region with >= 300 months of pre-1965 land precip
; boxtemp=boxes in each big region with >= 300 months of pre-1965 land air tem
; boxlists=boxes in each region with co-located MXD sites in them
;***********
;
;----------------------------------------------------
;
; Get chronology locations
;
restore,filename='allmxd.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
; Get region info
;
restore,'reg_mxdlists.idlsave'
;
; Get IPCC instrumental dataset grid
;
g=def_grid(/ipcc)
;
; For each region, locate the grid-boxes with data in
;
boxlists=fltarr(g.nx,g.ny,nreg)     ; boxes with MXD sites in
boxnumbs=fltarr(g.nx,g.ny,nreg)     ; number of MXD sites per box
print,'Locating boxes with MXD sites in'
for ireg = 0 , nreg-1 do begin
  print,regname(ireg),format='($,A20)'
  statx=statlon(treelist(0:ntree(ireg)-1,ireg))
  staty=statlat(treelist(0:ntree(ireg)-1,ireg))
  statval=fltarr(ntree(ireg))+1.
  fdout=gridit(g.nx,g.ny,g.x,g.y,statx,staty,statval,nstat=nstat)
  boxlists(*,*,ireg)=fdout(*,*)
  boxnumbs(*,*,ireg)=nstat(*,*)
endfor
;
; Get precipitation field (1961-90 means contain the land mask, since
; we're using the New and Hulme dataset which is complete by interpolation)
; *** IN FACT, IT DOES HAVE SOME MISSING VALUES IN IT, WHICH HAVE BEEN
; SET TO THE 61-90 MEAN.  I HAVE NOW RESET THESE TO MISSING, SO I HAVE
; AN INCOMPLETE DATASET ***
;
restore,'/cru/u2/f055/data/obs/grid/surface/precip_new_19611990.idlsave'
;  g,ltm5,nmon
;
; First identify the grid boxes in each region (since the Jones land
; air temperature is on the same grid, this will apply to that too)
;
boxregs=fltarr(g.nx,g.ny,nreg)     ; boxes in each region
x=fltarr(g.nx,g.ny)
for i = 0 , g.ny-1 do x(*,i)=g.x(*)
y=fltarr(g.nx,g.ny)
for i = 0 , g.nx-1 do y(i,*)=g.y(*)
print,'Locating boxes in each region'
for ireg = 0 , nreg-1 do begin
  print,regname(ireg)
  case ireg of
    0: begin         ; Northern Europe (NB shift east edge to get extra box!!!)
      kl=where( (x gt -15.) and (x le 63.) and (y gt 53.) and $
                ((x le 30.) or (y gt 60.)) , nkeep)
      end
    1: begin         ; Southern Europe
      kl=where( (x gt -15.) and (x le 30.) and (y le 53.) and (y gt 35.),nkeep)
      end
    2: begin         ; Northern Siberia
      kl=where( (x gt 61.) and (x le 128.) and (y gt 60.) , nkeep)
      end
    3: begin         ; Eastern Siberia (NB shift west edge to get extra box!!!)
      kl=where( ((x gt 127.) or (x le -170.)) and (y gt 60.) , nkeep)
      end
    4: begin         ; Central Asia
      kl=where( (x gt 30.) and (x le 115.) and (y le 60.) and (y gt 40.),nkeep)
      end
    5: begin         ; Tibetan Plateau
      kl=where( (x gt 60.) and (x le 115.) and (y le 40.) and (y gt 23.),nkeep)
      end
    6: begin         ; Western North America
      kl=where( (x gt -135.) and (x le -100.) and (y le 50.) and (y gt 30.) , $
                nkeep)
      end
    7: begin         ; North western North America
      kl=where( (x gt -170.) and (x le -105.) and (y gt 50.) and $
        ((x le -120.) or (y le 55.)) and ((x le -125.) or (y le 60.)) , nkeep)
      end
    8: begin         ; Eastern & central Canada
      kl=where( (x gt -125.) and (x le -45.) and (y gt 40.) and $
                ((x gt -120.) or (y gt 60.)) and $
                ((x gt -105.) or (y gt 55.)) and $
                ((x gt -100.) or (y gt 50.)) , nkeep)
      end
  endcase
  onefd=fltarr(g.nx,g.ny)*!values.f_nan
  onefd(kl)=1.
  boxregs(*,*,ireg)=onefd(*,*)
endfor
;
; Next identify those boxes in the region that have at least 300 months
; of pre-1965 land precipitation data in them
;
print,'Locating boxes with precipitation data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/precip_new_19011995.mon.nc')
tsmon=crunc_rddata(ncid,tst=[1901,0],tend=[1964,11],grid=g)
ncdf_close,ncid
tscount=total(finite(tsmon),3)
loadct,39
qtv,reverse(tscount,2),/scale
kl=where(tscount ge 300,nkeep)
print,'Globally, boxes with >= 300 months of pre-1965 data =',nkeep
oneprec=fltarr(g.nx,g.ny)*!values.f_nan
oneprec(kl)=1.
boxprec=fltarr(g.nx,g.ny,nreg)
for ireg = 0 , nreg-1 do boxprec(*,*,ireg)=boxregs(*,*,ireg)*oneprec(*,*)
;
; Now identify those boxes in the region that have at least 300 months
; of pre-1965 land air temperature data in them
;
print,'Locating boxes with land temperature data'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_jones_18511998.mon.nc')
tsmon=crunc_rddata(ncid,tst=[1851,0],tend=[1964,11],grid=g)
ncdf_close,ncid
tscount=total(finite(tsmon),3)
kl=where(tscount ge 300,nkeep)
print,'Globally, boxes with >= 300 months of pre-1965 data =',nkeep
onetemp=fltarr(g.nx,g.ny)*!values.f_nan
onetemp(kl)=1.
boxtemp=fltarr(g.nx,g.ny,nreg)
for ireg = 0 , nreg-1 do boxtemp(*,*,ireg)=boxregs(*,*,ireg)*onetemp(*,*)
;
save,filename='reg_mxdboxes.idlsave',nreg,boxlists,boxnumbs,regname,$
  g,boxregs,boxprec,boxtemp
;
end