; read in each of the files
datad=dd+'/data'
pc_read_var,object=ff,datadir=datad
pc_read_dim,object=dim,datadir=datad
pc_read_grid,object=grid,datadir=datad
z=grid.z
l1=dim.l1
m1=dim.m1
n1=dim.n1
l2=dim.l2
m2=dim.m2
n2=dim.n2
nz=dim.nz
nghostz=dim.nghostz
mean_rho=0.*findgen(nz)
mean_nd=0.*findgen(nz)
mean_cc=0.*findgen(nz)
; fix parameters ----
pc_read_ts,object=ts,datadir=datad
urms=mean(ts.urms[it1:it2])
Re=urms/(nu*kf)
;---------------------
for n=0,nz-1 do begin
  mean_rho[n]=mean(ff.lnrho[l1:l2,m1:m2,n+nghostz])
  mean_nd[n]=mean(ff.nd[l1:l2,m1:m2,n+nghostz])
  mean_cc[n]=mean(ff.cc[l1:l2,m1:m2,n+nghostz])
end
total_density=mean(mean_rho)
total_sdensity=mean(mean_cc)
total_pdensity=mean(mean_nd)
plot,mean_cc/total_sdensity,z[n1:n2],$
   xtitle='density',ytitle='z',title=dd,$
   xstyle=1,ystyle=1,xr=[0.8,2.5],ps=-4
xaxis=0.8+((2.5-0.8)/20.)*findgen(50)
zero=0.*findgen(50)
oplot,xaxis,zero
oplot,mean_nd/total_pdensity,z[n1:n2],ps=-2
oplot,mean_rho/total_density,z[n1:n2],ps=-1
legend,['fluid density','scalar density','dust density'],psym=[1,4,2],/right,$
      charsize=1.2
com=dd+"=create_struct('z',z,'n1',n1,'n2',n2,'urms',urms,'Re',Re,'nu',nu,'nud',nud,'betad',betad,'kf',kf,'rhomz',mean_rho,'rhopmz',mean_cc,'rhodmz',mean_nd,'rhom',total_density,'rhop',total_sdensity,'rhod',total_pdensity)"
void=execute(com)