;
restore,filename='reglists.idlsave'
;
loadct,0
def_1color,r,g,b,10,color='mlgrey'
multi_plot,nrow=4
if !d.name eq 'X' then window,ysize=850
;
; Ignore 1991 and 1992 (last two values) from all timeseries
;
bs=0.4
for i = 0 , nreg-1 do begin
  ;
  pause
  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 ge 1881,nkeep)
  y=wdens(keeplist)
  x=x(keeplist)
  yhist=tim_hist(y,binsize=bs,binvals=binvals)
  binvals=[binvals(0)-bs,binvals]  &  yhist=[0.,yhist]
  cpl_barts,binvals,yhist,$
    bar_color=10.,$
    title='Distribution of normalised anomalies for region: '+$
      regname(i),$
    ystyle=1,yrange=[0.,0.28],xrange=[-10,10],$
    xtitle='Density anomaly',$
    ytitle='Frequency'
  ;
  restore,filename='densadj_'+regname(i)+'.idlsave'
  mknormal,densadj,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
  keeplist=where(x lt 1993,nkeep)
  y=densadj(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
  ;
  restore,filename='rwidadj_'+regname(i)+'.idlsave'
  wdens=rwidadj
  mknormal,wdens,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
  keeplist=where(x lt 1993,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=1
  ;
  legend,['Instrumental temperature','Density','Ring width'],$
    colors=[10,!p.color,!p.color],thick=[10,5,1]
  ;
endfor
;
end