;
; Creates a regional timeseries and adjust it for changing sample size
; JUST DOES ALL, AND DOES TIME-DEPENDENT RBAR!!!
;
plot,[0,1]
multi_plot,nrow=4,layout='large'
if !d.name eq 'X' then begin
  window,ysize=900
  initx
endif
;
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='rbar_dens_18811960.idlsave'
restore,filename='reglists.idlsave'
;
i=0  ; do ALL only
;for i = 0 , nreg-1 do begin
;
  dens=density(*,treelist(0:ntree(i)-1,i))
  wt=weight(*,treelist(0:ntree(i)-1,i))
;
; Make a time-dependent rbar first
;
nt=n_elements(x)
rbar_td=fltarr(nt)*!values.f_nan
print,nt
for it = 0 , nt-1 do begin
  mask = where( finite(dens(it,*)) , nsamp)
  print,it,nsamp
  if nsamp gt 0 then begin
    if nsamp eq 1 then begin
      rbar_td(it)=1.               ; if n=1 then n'=1 regardless of rbar
    endif else begin
      ; get reduced grandr
      redr=reform(grandr(mask,*))
      redr=redr(*,mask)
      ; Total it and average it, but only do i<>j.
      totr=0.
      totn=0.
      for j = 1 , nsamp-1 do begin
        for k = 0 , j-1 do begin
          totr=totr+redr(j,k)        ; should have no missing values I hope
          totn=totn+1.
          if redr(j,k) eq 1. then message,'Ooops!'
        endfor
      endfor
      rbar_td(it)=totr/totn
    endelse
  endif
endfor
  ;
  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)
xadj=mkeffective(wdens,n,rbar=rbar_td,neff=neff)
wdens=xadj
;
; Now normalise w.r.t. 1881-1960
;
  mknormal,wdens,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
print,refmean,refsd
  ;
  ; Now plot them
  ;
  pause
  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
  ;
  plot,x,n,ytitle='Sample size',/xstyle
  ;
  plot,x,rbar_td,/xstyle,ytitle='rbar'
  plot,x,neff,/xstyle,ytitle='Effective sample size'
  ;
  ; Now save the weighted mean timeseries, for further analysis
  ;
  densadj=wdens
  save,filename='densadj_'+regname(i)+'.idlsave',x,densadj,n,neff,ncore
  ;
;endfor
;
end