;
; Need to provide a list of years in vector yrlist
;
yrlist=[1453,1601,1641,1695,1816]
md=[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=3,ncol=2,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,$
    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
;
; First plot a map of the regions and the chronology locations (with different
; colours according to chronology length
;
inter_boxfd,/nodata,map=map,labels=labels
;
!p.ticklen=0.03
lon_polar,map=map,[-180,-150.,-120.,-90.,-60.,-30.,30.,60.,120.,150.],/dotick
;
x1=[180.,60.,-40.,-105.,-105.,-107.]
y1=[30., 30., 30.,  30.,  55.,  55.]
x2=[180.,60.,-40.,-105.,-107.,-107.]
y2=[90., 90., 90.,  55.,  55.,  90.]
n=n_elements(x1)
for i = 0 , n-1 do begin
  map_plots,[x1(i),x2(i)],[y1(i),y2(i)],thick=6
endfor
;
xs=findgen(99)-40.
ys=replicate(54.8,99)
map_plots,xs,ys,thick=6
;
; Plots maps of chronology location
;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_varget,ncid,'latitude',ylat
ncdf_varget,ncid,'longitude',xlon
ncdf_varget,ncid,'year',timey
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_close,ncid
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
restore,filename='reglists.idlsave'
iall=where(regname eq 'ALL')
i=iall(0)
mxdall=density(*,treelist(0:ntree(i)-1,i))
x=xlon(treelist(0:ntree(i)-1,i))
y=ylat(treelist(0:ntree(i)-1,i))
;
allcut=reverse([1500,1600,1700,1800,1910])
ncut=n_elements(allcut)
;
cpl_usersym,/circle,/fill
for icut = 0 , ncut-1 do begin
  ct=allcut(icut)
  yw=where(timey eq ct)
  yw=yw(0)
  kl=where(finite(mxdall(yw,*)),nkeep)
  print,icut,ct,nkeep,90+icut
  if ct eq 1500 then print,x(kl)
  plots,x(kl),y(kl),psym=8,color=90+icut,symsize=!p.charsize*0.9
endfor
;
; 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
  densfd=rd1yr_norm(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,'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,-150.,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!Canomaly'
device,font_size=9
scale_vert,/off
!p.region=[0,0,0,0]
;
end