;
; Reads in instrumental 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
;
restore,filename='nc.instrboxes.idlsave'
thalf=10.
nstat=n_elements(boxlon)
;
fns='instr1978.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)
print,'nstat=',nstat
for i = 0 , nstat-1 do begin
  ts1=reform(boxts(i,*,*))              ; to give a mon , yr array
  totval=total(ts1(3:8,*),1,/nan)
  totn=total(finite(ts1(3:8,*)),1)
  dummy=where((totn gt 0) and (x ge 1901.) and (x le 1940.),nkkkk)
  print,i,nkkkk
  if nkkkk gt 10 then begin
    keeplist=where(totn gt 0,nkeep)
    ts2=totval*!values.f_nan
    ts2(keeplist)=totval(keeplist)/float(totn(keeplist))
    mknormal,ts2,x,refperiod=[1901,1940]
    filter_cru,thalf,tsin=ts2,tslow=ts3,/nan
    statvals(i)=ts3(iwant)
  endif else begin
    statvals(i)=!values.f_nan
  endelse
endfor
save,filename=fns,statvals
endelse
;
; Remove missing data and data south of 30N
;
keeplist=where(finite(statvals) and (boxlat ge 30.),nkeep)
if nkeep gt 0 then begin
   statvals=statvals(keeplist)
   boxlat=boxlat(keeplist)
   boxlon=boxlon(keeplist)
;  frac=frac(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=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)
;
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
;
;inter_const,gridfd,gx,gy,map=map,$
scale_vert,levels=levels,c_colors=levcol
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) Apr-Sep normalised temperature anomaly'
contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot
scale_vert,/off
;
end