;
; Reads in density data, puts a low-frequency filter through each retained
; grid-box, and plots anomaly (normalised w.r.t. 1901-40) map
;
; Get display ready
;
loadct,39
multi_plot,nrow=1,layout='large'
if !d.name eq 'X' then begin
  window,ysize=700,xsize=700
  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
;sm.thresh=0.1
;
; Read in data
;
ncid=ncdf_open('tree_rwid_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'width',density
ncdf_attget,ncid,'width','missing_value',valmiss
ncdf_varget,ncid,'fraction',frac
ncdf_varget,ncid,'longitude',statlon
ncdf_varget,ncid,'latitude',statlat
ncdf_close,ncid
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
thalf=10.
;
fns='rwid1978.idlsave'
dummy=findfile(fns,count=nfile)
if nfile gt 0 then begin
  restore,filename=fns
endif else begin
;
; For each available box, normalise w.r.t. 1901-40, and filter with thalf
;
statvals=fltarr(nstat)
iwant=where(x eq 1978)
fracvals=frac(iwant,*)
print,'nstat=',nstat
for i = 0 , nstat-1 do begin
  ts1=reform(density(*,i))
  dummy=where(finite(ts1) and (x ge 1901.) and (x le 1940.),nkkkk)
  print,i,nkkkk
  if nkkkk gt 10 then begin
    mknormal,ts1,x,refperiod=[1901,1940]
    filter_cru,thalf,tsin=ts1,tslow=ts3,/nan
    statvals(i)=ts3(iwant)
  endif else begin
    statvals(i)=!values.f_nan
  endelse
endfor
save,filename=fns,statvals,fracvals
endelse
;
; Remove missing data and data south of 30N
;
keeplist=where(finite(statvals) and (statlat ge 30.),nkeep)
if nkeep gt 0 then begin
   statvals=statvals(keeplist)
   statlat=statlat(keeplist)
   statlon=statlon(keeplist)
;   boxlat=boxlat(keeplist)
;   boxlon=boxlon(keeplist)
  fracvals=fracvals(keeplist)
;  densfd=densfd(keeplist)
;  x=x(keeplist)
;  y=y(keeplist)
endif
print,'No. of boxes 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=2.5  &  dy=2.5
gnx=360./dx
gny=(90.-30.)/dy
gx=findgen(gnx)*dx-180.
gy=90.-findgen(gny)*dy
gridfd=gridit(gnx,gny,gx,gy,statlon,statlat,statvals,fracvals,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)
;
levels=findgen(16)*0.4-3.
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
def_1color,r,g,b,10,color='purple'
def_1color,r,g,b,16,color='lgreen'
def_smearcolor,r,g,b,fromto=[10,16]
def_1color,r,g,b,17,color='lgrey'
def_1color,r,g,b,18,color='lyellow'
def_1color,r,g,b,24,color='red'
def_smearcolor,r,g,b,fromto=[18,24]
levcol=indgen(15)+10
;
scale_vert,levels=levels,c_colors=levcol
;
inter_const,gridfd,gx,gy,map=map,$
;inter_const,statvals,boxlon,boxlat,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,$
  miss_grey='white',$
  levels=levels,c_colors=levcol,$
  /cell_fill,fdsm=fdsm,$
  title='1978 (20 yr) normalised width anomaly'
;
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot
;
scale_vert,/off
;
end