;
; Need to provide a list of years in vector yrlist
; In fact, yrlist is 2D array (n,m), with n different maps each made up
; of the mean of m different years.  If m equal NaN then it is ignored.
; titlist=list of plot titles
;
if n_elements(yrlist) eq 0 then message,'Need list of years in yrlist'
;
; Get display ready
;
loadct,39
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
;
lsz=size(yrlist)
nmap=lsz(2)
for i = 0 , nmap-1 do begin
  alist=reform(yrlist(i,*))
  onelist=alist( where( finite(alist) ) )
  nyr=n_elements(onelist)
  for iyy = 0 , nyr-1 do begin
    iyr=onelist(iyy)
    print,iyy,iyr
    onefd=rd1yr(iyr,x=x,y=y,frac=frac)
    if iyy eq 0 then begin
      nchron=n_elements(onefd)
      allfd=fltarr(nchron,nyr)
    endif
    allfd(*,iyy)=onefd
  endfor
  misslist=where(allfd eq -9.99,nmiss)
  if nmiss gt 0 then allfd(misslist)=!values.f_nan
  summap=total(allfd,/nan,2)
  coumap=total(finite(allfd),2)
  densfd=summap/float(coumap)
;
; Remove missing data and data south of 30N
;
keeplist=where(finite(densfd) 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,'No. of points with data:',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.
levels=levels*0.5
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=titlist(i)
;
endfor
;
end