;
; Creates a regional timeseries and adjust it for changing sample size
;
plot,[0,1]
multi_plot,nrow=4
if !d.name eq 'X' then begin
  window,ysize=900
  initx
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
;
; Get regional tree lists and rbar
;
restore,filename='reglists.idlsave'
restore,filename='rbar_dens_full.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 need to adjust for varying sample size
;
xadj=mkeffective(wdens,n,rbar=allrbar(i),neff=neff)
wdens=xadj
;
; Now normalise w.r.t. 1901-1940
;
help,wdens
  mknormal,wdens,x,refperiod=[1901,1940],refmean=refmean,refsd=refsd
help,wdens
print,refmean,refsd
  ;
  ; Now plot them
  ;
  pause
;  plot,x,wdens,psym=10,title='Weighted mean normalised ring width anomaly'+$
;    ' for region: '+regname(i),xtitle='Year'
  filter_cru,20,tsin=wdens,tslow=tslow,/nan 
  cpl_barts,x,wdens,title='Weighted mean normalised max density anomaly'+$ 
    ' for region: '+regname(i),xtitle='Year',/xstyle,$ 
    zeroline=tslow,yrange=[-8,4]
  oplot,x,tslow,thick=2 
  moms=moment(wdens(where(finite(wdens))),sdev=sdev) 
  oplot,!x.crange,[0.,0.]
  oplot,!x.crange,[-2.,-2.],linestyle=1
;  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 
;  plot,x,neff,ytitle='Effective sample size',/xstyle
  ;
  ; Now save the weighted mean timeseries, for further analysis
  ;
  densadj=wdens
  save,filename='densadj_'+regname(i)+'.idlsave',x,densadj,n,neff,ncore
  ;
endfor
;
end