;
; Get display ready
;
loadct,39
multiplot,1
;
; 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.
;
; Small amount of smoothing, with filling in around the edges to extend to
; region that is plotted
;
sm=def_sm()
sm.thresh=0.1
;
print,'Enter year required:'
read,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
;
inter_const,gridfd,gx,gy,map=map,$
maxdist=3000.,$
gs=[2.,2.],$
sm=sm,labels=labels,$
shade=1,sh_thresh=0.,$
levels=findgen(21)*0.5-5.,/follow,$
title=string(iyr,format='(I4)'),$
xtitle='Tree ring density anomaly'
;
end