;--------Here the unit are SI unit-----------

;grid=128
;grid_len=1.0/grid
grid_len=L/grid
k_prefactor=3.13
density_ratio=1000.0
g=9.80
good_t=50 ;seconds
;read initial parameters------
pc_read_param,obj=start
boxsize=abs(start.xyz0*2)
boxsize=boxsize(0)
a0=start.ap0

pc_read_param,/param2,obj=param
nu=param.nu*1.0
k1=2*!pi/boxsize
kf=k1*k_prefactor

pc_read_ts,obj=ts
good=where(ts.t GT good_t)
good=good[0]
urms=mean(ts.urms(good:*))
epsk=mean(ts.epsk(good:*))

print,'boxsize=',boxsize
print,'urms=',urms
print,'epsk=',epsk
print,'nu=',nu
print,'k1=',k1
print,'kf=',kf
;Komogorov scales
;Re=urms/(nu*kf)
Re=urms*boxsize/(nu*k_prefactor) ; Re=Urms*L_inj/nu
eta=boxsize/Re^0.75
scale_seperation=(boxsize/k_prefactor)/eta ;note that the largest scale in turbulence is the energy injection scale
tau=(boxsize/urms)/sqrt(Re)
epson=urms^2/tau
v_eta=(nu*epson)^0.25
;t_turnover=1.0/(kf*urms)  ; T_eddy=L_inj/U_rms=L/(k_prefactor*U_rms)
t_turnover=boxsize/(k_prefactor*urms)
nu_tur=urms/(3.0*kf)
t_tur=boxsize^2/nu_tur
Re_mesh=urms*grid_len/nu     ;Re_mesh shoud be <=5.
print,'Re=',Re
print,'eta=',eta
print,'tau=',tau
print,'epson=',epson
print,'v_eta=',v_eta
print,'t_turnover=',t_turnover
print,'nu_tur=',nu_tur
print,'t_turDiffusion=',t_tur
print,'Re_mesh=',Re_mesh
print,'scale_seperation',scale_seperation
;stokes time
tau_p=2.0/9.0*density_ratio*a0^2/nu
v_p=tau_p*g
St=tau_p/tau
v_ratio=v_p/v_eta

print,'tau_p=',tau_p
print,'v_p=',v_p
print,'St=',St
print,'v_ratio=',v_ratio
END