;
; Need to provide a list of years in vector yrlist
;
yrlist=[1783,1816]
if n_elements(yrlist) eq 0 then message,'Need list of years in yrlist'
;
; Get display ready
;
!y.omargin=[0,1]
loadct,39
def_1color,r,g,b,40,color='purple'
def_1color,r,g,b,42,color='dblue'
def_1color,r,g,b,43,color='blue'
def_1color,r,g,b,46,color='lgreen'
def_1color,r,g,b,47,color='lyellow'
def_1color,r,g,b,50,color='red'
def_smearcolor,r,g,b,fromto=[40,42]
def_smearcolor,r,g,b,fromto=[43,46]
def_smearcolor,r,g,b,fromto=[47,50]
multi_plot,nrow=2,ncol=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=800,xsize=800
  fac=1.
endif else begin
  fac=2.
endelse
;
; N. Hemi. map down to 30N, no labels
;
map=def_map()
map.limit=[0,-180,90,180] & map.origx=0.
labels=def_labels(/off)
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(/off)
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]
nyr=n_elements(yrlist)
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=5.  &  dy=5.
gnx=360./dx
gny=(90.-0.)/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=[-12,-6,-4,-3,-2,-1,-0.5,0,0.5,1,2,3]
nlev=n_elements(levels)
zeroloc=where(levels eq 0,nzero)
c_labels=intarr(nlev)+0  &  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.
c_colors=findgen(11)+40.
print,levels
print,c_thick
;
print,min(gridfd,/nan,max=mmm) & print,mmm
inter_boxfd,gridfd,gx,gy,map=map,$
  sm=sm,labels=labels,fdsm=fdsm,$
  miss_grey='white',$
  levels=levels,c_colors=c_colors,$
  title=string(iyr,format='(I4)'),/scale
;  xtitle='Tree ring density anomaly'
;
endfor
;
end