;
restore,filename='reglists.idlsave'
;
multi_plot,nrow=4
;
; Ignore 1991 and 1992 (last two values) from all timeseries
;
bs=0.4
for i = 0 , nreg-1 do begin
  restore,filename='densadj_'+regname(i)+'.idlsave'
  mknormal,densadj,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
  keeplist=where(x lt 1991,nkeep)
  y=densadj(keeplist)
  x=x(keeplist)
  yhist=tim_hist(y,binsize=bs,binvals=binvals)
  pause
  binvals=[binvals(0)-bs,binvals]  &  yhist=[0.,yhist]
  plot,binvals,yhist,psym=10,$
    title='Distribution of normalised anomalies for region: '+$
      regname(i),$
    ystyle=1,yrange=[0.,0.28],xrange=[-10,10],$
    xtitle='Density anomaly',$
    ytitle='Frequency'
  legend,['Instrumental temperature','Density','Ring width'],$
    linestyle=[0,0,2],thick=[5,1,3]
  restore,filename='rwidadj_'+regname(i)+'.idlsave'
  wdens=rwidadj
  mknormal,wdens,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
  keeplist=where(x lt 1991,nkeep)
  y=wdens(keeplist)
  x=x(keeplist)
  yhist=tim_hist(y,binsize=bs,binvals=binvals)
  binvals=[binvals(0)-bs,binvals]  &  yhist=[0.,yhist]
  oplot,binvals,yhist,psym=10,linestyle=2,thick=3
  restore,filename='instrts_'+regname(i)+'.idlsave'
  totval=total(winstr(*,3:8),2,/nan) 
  totn=total(finite(winstr(*,3:8)),2) 
  keeplist=where(totn gt 0,nkeep) 
  wdens=totval*!values.f_nan                       
  wdens(keeplist)=totval(keeplist)/float(totn(keeplist)) 
  mknormal,wdens,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
  keeplist=where(x lt 1991,nkeep)
  y=wdens(keeplist)
  x=x(keeplist)
  yhist=tim_hist(y,binsize=bs,binvals=binvals)
  binvals=[binvals(0)-bs,binvals]  &  yhist=[0.,yhist]
  oplot,binvals,yhist,psym=10,thick=5
endfor
;
end