;
loadct,39
multiplot,1
if !d.name eq 'X' then begin
  window,xsize=950,ysize=950
  chf=0.75
endif else begin
  device,xoffset=1,xsize=19,yoffset=2,ysize=25
  chf=0.6
endelse
;
map=def_map(/npolar)
map.limit(0)=30.
labels=def_labels(/off)
;
print,'Enter year required:'
read,iyr
densfd=rd1yr(iyr,x=x,y=y,frac=frac)
;
; Remove missing data and data south of 30N
;
keeplist=where((densfd ne -9.99) and (y ge 30.),nkeep)
if nkeep gt 0 then begin
  frac=frac(keeplist)
  densfd=densfd(keeplist)
  x=x(keeplist)
  y=y(keeplist)
endif
;
; First of all, store each station value into its 0.5 by 0.5 grid box.
; When more than one falls in a box, average them - this is where the
; weighting by the fraction of cores available comes into it!  It also
; prevents duplicate points going forward, which can upset spherical
; triangulation.
;
dx=1.  &  dy=1.
gnx=360./dx
gny=(90.-30.)/dy
gx=findgen(gnx)*dx-180.
gy=90.-findgen(gny)*dy
gridfd=gridit(gnx,gny,gx,gy,x,y,densfd,frac,nstat=chron,minv=minv,maxv=maxv)
;
inter_boxfd,gridfd,gx,gy,map=map,$
  /nodata,$
  labels=labels,$
  title=string(iyr,format='(I4)'),$
  xtitle='Weighted mean tree ring density anomaly * 10'
;
for i = 0 , gnx-1 do begin
  for j = 0 , gny-1 do begin
    if finite(gridfd(i,j)) then $
      xyouts,gx(i),gy(j),string(gridfd(i,j)*10.,format='(I3)'),$
        alignment=0.5,charthick=0.5,charsize=chf*!p.charsize
  endfor
endfor
;
pause
;
inter_boxfd,gridfd,gx,gy,map=map,$
  /nodata,$
  labels=labels,$
  title=string(iyr,format='(I4)'),$
  xtitle='Maximum tree ring density anomaly * 10'
;
for i = 0 , gnx-1 do begin
  for j = 0 , gny-1 do begin
    if finite(maxv(i,j)) then $
      xyouts,gx(i),gy(j),string(maxv(i,j)*10.,format='(I3)'),$
        alignment=0.5,charthick=0.5,charsize=chf*!p.charsize
  endfor
endfor
;
pause
;
inter_boxfd,gridfd,gx,gy,map=map,$
  /nodata,$
  labels=labels,$
  title=string(iyr,format='(I4)'),$
  xtitle='Minimum mean tree ring density anomaly * 10'
;
for i = 0 , gnx-1 do begin
  for j = 0 , gny-1 do begin
    if finite(minv(i,j)) then $
      xyouts,gx(i),gy(j),string(minv(i,j)*10.,format='(I3)'),$
        alignment=0.5,charthick=0.5,charsize=chf*!p.charsize
  endfor
endfor
;
end