;
; Need to provide a list of years in vector yrlist
;
nyr=n_elements(yrlist)
if nyr eq 0 then message,'Need list of years in yrlist'
;
; Get display ready
;
loadct,39
if nyr eq 2 then multi_plot,nrow=2,ncol=1 $
            else multi_plot,nrow=3,ncol=2,layout='large'
if !d.name eq 'X' then begin
  window,ysize=950
  fac=1.
endif else begin
  fac=2.
endelse
;
; 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.
coast=def_coast(/get_device)
coast.double=0
;
; Small amount of smoothing, with filling in around the edges to extend to
; region that is plotted
;
sm=def_sm()
sm.method=2
sm.thresh=1.0
;
;print,'Enter year required:'
;read,iyr
;yrlist=[1601,1641,1669,1699,1816,1817,1818,1878,1884,1912,1964]
;yrlist=[1818,1878,1884,1912,1964]
for iyy = 0 , nyr-1 do begin
  iyr=yrlist(iyy)
  print,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
;
levels=findgen(21)-10.
nlev=n_elements(levels)
zeroloc=where(levels eq 0,nzero)
c_labels=intarr(nlev)+1  &  if nzero gt 0 then c_labels(zeroloc)=0
c_charsize=0.6
;c_thick=(levels lt 0.)*fac+2.  &  if nzero gt 0 then c_thick(zeroloc)=1.
c_thick=1.5*fac  &  if nzero gt 0 then c_thick(zeroloc)=1.
print,levels
print,c_thick
;
inter_const,gridfd,gx,gy,map=map,$
  maxdist=3000.,$
  gs=[2.,2.],$
  sm=sm,labels=labels,$
  shade=1,sh_thresh=0.,sh_grey=[235,180],miss_grey='white',$
  levels=levels,c_labels=c_labels,c_charsize=c_charsize,c_thick=c_thick,$
  /follow,$
  title=string(iyr,format='(I4)')
;  xtitle='Tree ring density anomaly'
;
endfor
;
end