program kolmogorov_scales
    implicit none
    double precision, parameter :: C_mu = 0.09d0, C_K = 1.5d0

    double precision :: eps_in, nu_in, rho, mu, k_tke, U_mean, L_dom
    double precision :: nu, eps, u_rms
    double precision :: eta, tau_eta, u_eta
    double precision :: lam_g, lam_iso, Re_lam
    double precision :: L_int, T_int, u_int
    double precision :: L_eta, L_lam, lam_eta
    double precision :: Re_L, Re_dom
    double precision :: N_grid, N_time, E_L, E_eta
    integer :: i, n_pts, ios
    double precision :: ec, etc, ttc, lmc, Lc, rlc, nc

    read(*,*,iostat=ios) eps_in;  if(ios/=0)then;write(*,*)'ERROR: Invalid epsilon.';stop;end if
    read(*,*,iostat=ios) nu_in;   if(ios/=0)then;write(*,*)'ERROR: Invalid nu.';stop;end if
    read(*,*,iostat=ios) rho;     if(ios/=0)then;write(*,*)'ERROR: Invalid rho.';stop;end if
    read(*,*,iostat=ios) mu;      if(ios/=0)then;write(*,*)'ERROR: Invalid mu.';stop;end if
    read(*,*,iostat=ios) k_tke;   if(ios/=0)then;write(*,*)'ERROR: Invalid k.';stop;end if
    read(*,*,iostat=ios) U_mean;  if(ios/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
    read(*,*,iostat=ios) L_dom;   if(ios/=0)then;write(*,*)'ERROR: Invalid L_domain.';stop;end if

    if(eps_in<=0d0)then;write(*,*)'ERROR: epsilon must be > 0.';stop;end if
    if(k_tke<=0d0)then;write(*,*)'ERROR: k must be > 0.';stop;end if

    if (nu_in > 0d0) then
        nu = nu_in
    else
        if(rho<=0d0.or.mu<=0d0)then;write(*,*)'ERROR: Need nu>0 or rho,mu>0.';stop;end if
        nu = mu / rho
    end if
    eps = eps_in

    ! Kolmogorov scales
    eta     = (nu**3 / eps)**0.25d0
    tau_eta = sqrt(nu / eps)
    u_eta   = (nu * eps)**0.25d0

    ! RMS velocity
    u_rms = sqrt(2d0 * k_tke / 3d0)

    ! Taylor microscale
    lam_g   = sqrt(10d0 * nu * k_tke / eps)
    lam_iso = sqrt(15d0 * nu * u_rms**2 / eps)
    Re_lam  = u_rms * lam_g / nu

    ! Integral scales
    L_int = k_tke**1.5d0 / eps
    T_int = k_tke / eps
    u_int = u_rms

    ! Scale ratios
    L_eta   = L_int / eta
    L_lam   = L_int / lam_g
    lam_eta = lam_g / eta

    ! Reynolds numbers
    Re_L   = u_rms * L_int / nu
    if (L_dom > 0d0 .and. U_mean > 0d0) then
        Re_dom = U_mean * L_dom / nu
    else
        Re_dom = 0d0
    end if

    ! DNS estimation
    N_grid = L_eta**3
    N_time = T_int / tau_eta

    ! Energy
    E_L  = u_rms**2
    E_eta = sqrt(nu * eps)

    ! Output
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   KOLMOGOROV SCALES CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,ES14.6,A)') '  epsilon (dissipation)   = ', eps, ' m2/s3'
    write(*,'(A,ES14.6,A)') '  Kinematic Visc (nu)     = ', nu, ' m2/s'
    write(*,'(A,ES14.6,A)') '  k (TKE)                 = ', k_tke, ' m2/s2'
    write(*,'(A,F12.4,A)')  '  u_rms = sqrt(2k/3)      = ', u_rms, ' m/s'
    write(*,'(A,F12.4,A)')  '  U_mean                  = ', U_mean, ' m/s'
    write(*,'(A,F12.4,A)')  '  L_domain                = ', L_dom, ' m'
    write(*,*)

    write(*,'(A)') '--- KOLMOGOROV SCALES (smallest) ----------------------------'
    write(*,'(A,ES14.6,A)') '  eta (length)            = ', eta, ' m'
    write(*,'(A,ES14.6,A)') '  tau_eta (time)          = ', tau_eta, ' s'
    write(*,'(A,ES14.6,A)') '  u_eta (velocity)        = ', u_eta, ' m/s'
    write(*,'(A,F12.4)')    '  Re_eta (= 1 by def.)    = ', u_eta * eta / nu
    write(*,*)

    write(*,'(A)') '--- TAYLOR MICROSCALE ---------------------------------------'
    write(*,'(A,ES14.6,A)') '  lambda (general)        = ', lam_g, ' m'
    write(*,'(A,ES14.6,A)') '  lambda (isotropic)      = ', lam_iso, ' m'
    write(*,'(A,F14.4)')    '  Re_lambda               = ', Re_lam
    write(*,*)

    write(*,'(A)') '--- INTEGRAL SCALES (largest) -------------------------------'
    write(*,'(A,ES14.6,A)') '  L_int = k^3/2 / eps     = ', L_int, ' m'
    write(*,'(A,ES14.6,A)') '  T_int = k / eps          = ', T_int, ' s'
    write(*,'(A,F12.6,A)')  '  u_int ~ u_rms           = ', u_int, ' m/s'
    write(*,*)

    write(*,'(A)') '--- SCALE RATIOS --------------------------------------------'
    write(*,'(A,ES14.6)')   '  L / eta                 = ', L_eta
    write(*,'(A,F14.4)')    '  L / lambda              = ', L_lam
    write(*,'(A,F14.4)')    '  lambda / eta             = ', lam_eta
    write(*,'(A,ES14.6)')   '  T_int / tau_eta          = ', T_int / tau_eta
    write(*,'(A,F14.4)')    '  u_int / u_eta            = ', u_int / u_eta
    write(*,*)

    write(*,'(A)') '--- REYNOLDS NUMBERS ----------------------------------------'
    write(*,'(A,ES14.6)')   '  Re_L = u_rms*L/nu       = ', Re_L
    write(*,'(A,F14.4)')    '  Re_lambda               = ', Re_lam
    write(*,'(A,ES14.6)')   '  Re_domain               = ', Re_dom
    write(*,*)

    write(*,'(A)') '--- DNS ESTIMATION ------------------------------------------'
    write(*,'(A,ES14.6)')   '  N_grid ~ (L/eta)^3      = ', N_grid
    write(*,'(A,F14.2)')    '  (L/eta)^(1/3) per dim   = ', L_eta
    write(*,'(A,ES14.6)')   '  N_timesteps ~ T/tau_eta  = ', T_int / tau_eta
    write(*,'(A,ES14.6)')   '  Total cost ~ N*Nt       = ', N_grid * T_int / tau_eta
    write(*,*)

    write(*,'(A)') '--- ENERGY CASCADE ------------------------------------------'
    write(*,'(A,ES14.6,A)') '  E_L ~ u_rms^2           = ', E_L, ' m2/s2'
    write(*,'(A,ES14.6,A)') '  E_eta ~ sqrt(nu*eps)    = ', E_eta, ' m2/s2'
    write(*,'(A,F12.4)')    '  Kolmogorov constant C_K  = ', C_K
    write(*,'(A)')          '  E(kappa) = C_K * eps^(2/3) * kappa^(-5/3)'
    write(*,*)

    ! Profile
    write(*,'(A)') '--- PROFILE vs DISSIPATION RATE -----------------------------'
    write(*,'(A)') '  epsilon      eta          tau_eta      lambda       L_int        Re_lam      (L/eta)'
    write(*,'(A)') '  ---------------------------------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        ec = eps / 100d0 * 10d0**(3d0 * dble(i-1) / dble(n_pts-1))
        etc = (nu**3 / ec)**0.25d0
        ttc = sqrt(nu / ec)
        lmc = sqrt(10d0 * nu * k_tke / ec)
        Lc  = k_tke**1.5d0 / ec
        rlc = u_rms * lmc / nu
        nc  = Lc / etc

        write(*,'(ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,F10.2,2X,ES12.4)') &
            ec, etc, ttc, lmc, Lc, rlc, nc
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') "  eta = (nu^3/eps)^(1/4)     (Kolmogorov length)"
    write(*,'(A)') "  tau_eta = (nu/eps)^(1/2)   (Kolmogorov time)"
    write(*,'(A)') "  u_eta = (nu*eps)^(1/4)     (Kolmogorov velocity)"
    write(*,'(A)') "  lambda = sqrt(10*nu*k/eps)  (Taylor microscale)"
    write(*,'(A)') "  L = k^(3/2) / eps           (integral scale)"
    write(*,'(A)') "  Re_lambda = u_rms*lambda/nu"
    write(*,'(A)') "  N_dns ~ (L/eta)^3 ~ Re_L^(9/4)"
    write(*,'(A)') '============================================================'

end program kolmogorov_scales
