;
; Need to provide a list of years in vector yrlist
;
; TO MAKE GRIDIT WORK, NEED FIRST TO DO:
;
; .run /cru/u2/f055/tree5/gridit.pro
;
yrlist=indgen(12)+1809.
md=[2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300]
;
; Get display ready
;
loadct,39
def_1color,r,g,b,90,color='black'
def_1color,r,g,b,91,color='purple'
def_1color,r,g,b,92,color='blue'
def_1color,r,g,b,93,color='green'
def_1color,r,g,b,94,color='red'
;
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=4,ncol=3,layout='large'
if !d.name eq 'X' then begin
  window, ysize=950
  !p.font=1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9,ysize=24.2
;    xsize=17.5,xoffset=1.5
endelse
!x.omargin=[0,15]
;
; N. Hemi. map down to 30N, no labels
;
map=def_map(/npolar)
map.limit(0)=30.
labels=def_labels(/off)
;
; 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=0.05
;
; Get MXD data
;
print,'Reading MXD data'
restore,filename='allmxd.idlsave'
;  nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$
;  mxd,fraction,timey,nyr
;
; Now plot the yearly maps (read from file containing mxd normalised
; w.r.t. 1881-1960).
;
nyr=n_elements(yrlist)
for iyy = 0 , nyr-1 do begin
  iyr=yrlist(iyy)
  print,'Year:',iyr
  kl=where(timey eq iyr)
  densfd=reform(mxd(kl,*))
  frac=reform(fraction(kl,*))
  x=statlon
  y=statlat
  ;
  ; Remove missing data
  ;
  keeplist=where((densfd ne -9.99),nkeep)
  if nkeep gt 0 then begin
    frac=frac(keeplist)
    densfd=densfd(keeplist)
    x=x(keeplist)
    y=y(keeplist)
  endif
  print,'Number of chronology values:',nkeep
  ;
  ; First of all, store each station value into its 1. by 1. 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=2.  &  dy=2.
  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)
  print,'Min/Max values:',min(gridfd,/nan,max=mmm) & print,mmm
  ;
  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_colors=findgen(11)+40.
  ;
  inter_const,gridfd,gx,gy,map=map,$
    maxdist=md(iyy),$
    gs=[2.,2.],$
    sm=sm,labels=labels,fdsm=fdsm,$
    miss_grey='white',$
    /cell_fill,levels=levels,c_colors=c_colors
  contour,fdsm.fd,fdsm.x,fdsm.y,/overplot,levels=levels,c_labels=c_labels,$
    c_charsize=c_charsize,/follow
  ;
  xyouts,-145.,27.,string(iyr,format='(I4)'),align=1.
  ;
endfor
;
; Now plot a colour scale
;
device,font_size=12
!p.multi(0)=1
!p.region=[0.9,0.35,1.,0.65]
!y.margin=map.ymargin
scale_vert,levels=levels,$
  c_colors=c_colors,/noextremes,$
  title='   Normalised!Cdensity!C anomaly'
device,font_size=9
scale_vert,/off
!p.region=[0,0,0,0]
;
end