;
; Plots mean timeseries of density anomaly, no. of contributing chronologies
; and fraction of cores available
;
plot,[0,1]
multiplot,4
loadct,39
def_1color,r,g,b,1,color='red'
def_1color,r,g,b,2,color='blue'
if !d.name eq 'X' then begin
  initx
  window,ysize=900
endif else begin
  device,xoffset=2,xsize=17
endelse

;
ncid=ncdf_open('tree_dens_nh.nc')
ncdf_diminq,ncid,'time',dummy,ntime
ncdf_diminq,ncid,'station',dummy,nstat
ncdf_varget,ncid,'year',x
ncdf_varget,ncid,'density',density
ncdf_attget,ncid,'density','missing_value',valmiss
ncdf_varget,ncid,'fraction',weight
ncdf_close,ncid
;
misslist=where(density eq valmiss,nmiss)
density(misslist)=!values.f_nan
;
restore,filename='reglists.idlsave'
;
for i = 0 , nreg-1 do begin
  dens=density(*,treelist(0:ntree(i)-1,i))
  wt=weight(*,treelist(0:ntree(i)-1,i))
  ;
  n=total(dens*0.+1.,2,/nan)       ; compute timeseries of number of chronols
  ncore=total(dens*0.+wt,2,/nan)   ; compute timeseries of summed core fraction
  totwdens=total(dens*wt,2,/nan)
  totwvals=total(float(finite(dens))*wt,2)
  wdens=totwdens/float(totwvals)               ; this is the weighted anomaly
;
; Now normalise w.r.t. 1901-1970
;
;help,wdens
;  mknormal,wdens,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
;help,wdens
;print,refmean,refsd
  ;
  ; Now plot them
  ;
  pause
;  plot,x,wdens,psym=10,title='Weighted mean normalised density anomaly'+$
;    ' for region: '+regname(i),xtitle='Year'
  filter_cru,20,tsin=wdens,tslow=tslow,/nan
  cpl_barts,x,wdens,title='Weighted mean normalised density anomaly'+$
    ' for region: '+regname(i),xtitle='Year',/xstyle,$
    zeroline=tslow,yrange=[-4,4],noclip=1
  oplot,x,tslow,thick=2
  moms=moment(wdens(where(finite(wdens))),sdev=sdev)
  oplot,!x.crange,[moms(0),moms(0)]
  oplot,!x.crange,[moms(0)+sdev,moms(0)+sdev],linestyle=1
  oplot,!x.crange,[moms(0)-sdev,moms(0)-sdev],linestyle=1
  ;
  ; Now save the weighted mean timeseries, for further analysis
  ;
  save,filename='densts_'+regname(i)+'.idlsave',x,wdens,n,ncore
  ;
endfor
;
end