;
; Get display ready
;
loadct,39
multiplot,1
;
; N. Hemi. map down to 30N, no labels
;
map=def_map(/npolar)
map.limit(0)=30.
labels=def_labels(/off)
labels.gridon=1  &  labels.dlon=90.
;
; Small amount of smoothing, with filling in around the edges to extend to
; region that is plotted
;
sm=def_sm()
sm.thresh=0.1
;
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
print,nkeep
;
; 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)
;
; Convert boxes back to a list of stations
;
gx2d=gx # (intarr(gny)+1.)
gy2d=transpose(gy # (intarr(gnx)+1.))
keeplist=where(finite(gridfd))
gridfd=gridfd(keeplist)
gx=gx2d(keeplist)
gy=gy2d(keeplist)
;
pause
;
inter_const,gridfd,gx,gy,map=map,$
  maxdist=3000.,$
  gs=[2.,2.],$
  sm=sm,labels=labels,$
  shade=1,sh_thresh=0.,$
  levels=findgen(21)*0.5-5.,/follow,$
  title=string(iyr,format='(I4)'),$
  xtitle='Tree ring density anomaly'
;
end