;
; Creates a regional timeseries and adjust it for changing sample size
;
plot,[0,1]
multi_plot,nrow=6
if !d.name eq 'X' then begin
  window,ysize=900
  initx
endif else begin
  device,xoffset=2,xsize=17
endelse
;
ncid=ncdf_open('longvolc.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
wt=weight(*,*)
;
nms=['Polar Urals','Tornetrask','Jasper']
for i = 0 , 2 do begin
  cpl_barts,x,density(*,i),title='longvolc record for '+nms(i),$
    xtitle='Year',/xstyle,$
    zeroline=0.,yrange=[-4,4]
  plot,x,wt(*,i),title='effective sample size for '+nms(i),$
    xtitle='Year',/xstyle
endfor
;
restore,filename='rbar_longvolc.idlsave'
;
dens=density(*,*)
wt=weight(*,*)
;
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)
totudens=total(dens,2,/nan)
totwvals=total(float(finite(dens))*wt,2)
totuvals=total(float(finite(dens)),2)
wdens=totwdens/float(totwvals)               ; this is the weighted anomaly
udens=totudens/float(totuvals)               ; this is the unweighted anomaly
;
pause
cpl_barts,x,wdens,title='longvolc record: simple weighted mean',$
  xtitle='Year',/xstyle,$ 
  zeroline=0.,yrange=[-4,4]
cpl_barts,x,udens,title='longvolc record: simple mean',$
  xtitle='Year',/xstyle,$ 
  zeroline=0.,yrange=[-4,4]
;
; Now need to adjust for varying sample size
;
xadj=mkeffective(udens,n,rbar=allrbar,neff=neff,/local)
keepold=udens
udens=xadj
;
; Now normalise w.r.t. 1901-1970
;
mknormal,keepold,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
mknormal,udens,x,refperiod=[1901,1970],refmean=refmean,refsd=refsd
print,refmean,refsd
;
; Now plot them
;
;filter_cru,20,tsin=keepold,tslow=tslow,/nan 
cpl_barts,x,keepold,title='Mean record normalised over 1901-70',$
  xtitle='Year',/xstyle,$ 
  zeroline=tslow,yrange=[-6,6]
;  oplot,x,tslow,thick=2 
oplot,!x.crange,[0.,0.]
oplot,!x.crange,[-2.,-2.],linestyle=1
oplot,!x.crange,[2.,2.],linestyle=1
;
;filter_cru,20,tsin=udens,tslow=tslow,/nan 
cpl_barts,x,udens,title='Mean record adjusted for sample size and then normalised over 1901-70',$
  xtitle='Year',/xstyle,$ 
  zeroline=0.,yrange=[-6,6]
;oplot,x,tslow,thick=2 
;moms=moment(udens(where(finite(udens))),sdev=sdev) 
oplot,!x.crange,[0.,0.]
oplot,!x.crange,[-2.,-2.],linestyle=1
oplot,!x.crange,[2.,2.],linestyle=1
;
plot,x,n,title='Number of chronologies',yrange=[-0.5,3.5],/xstyle
plot,x,neff,title='Effective number of independent chronologies',$
  yrange=[-0.5,3.5],/xstyle
;
; Now save the weighted mean timeseries, for further analysis
;
densadj=udens
save,filename='longvolcadj.idlsave',x,densadj,n,neff,ncore
;
end