program turbulent_ke
    implicit none
    double precision, parameter :: C_mu = 0.09d0

    double precision :: up, vp, wp, U_mean, delta, rho, mu, eps_in, L_in
    double precision :: nu, k_val, rho_k, u_rms
    double precision :: TI, TI_pct
    double precision :: eps_val, L_val, lambda_val
    double precision :: eta, tau_eta, u_eta
    double precision :: nu_t_ke, nu_t_kw, omega_val, nut_ratio
    double precision :: Re_t, Re_lam
    double precision :: b11, b22, b33, II_inv, III_inv
    double precision :: P_est
    logical :: isotropic_flag
    integer :: i, n_pts, iostat_val
    double precision :: ti_cur, k_c, eps_c, nut_c, ret_c, eta_c, urms_c

    ! ---- read --------------------------------------------------------
    read(*,*,iostat=iostat_val) up;      if(iostat_val/=0)then;write(*,*)'ERROR: Invalid u_prime.';stop;end if
    read(*,*,iostat=iostat_val) vp;      if(iostat_val/=0)then;write(*,*)'ERROR: Invalid v_prime.';stop;end if
    read(*,*,iostat=iostat_val) wp;      if(iostat_val/=0)then;write(*,*)'ERROR: Invalid w_prime.';stop;end if
    read(*,*,iostat=iostat_val) U_mean;  if(iostat_val/=0)then;write(*,*)'ERROR: Invalid U_mean.';stop;end if
    read(*,*,iostat=iostat_val) delta;   if(iostat_val/=0)then;write(*,*)'ERROR: Invalid delta.';stop;end if
    read(*,*,iostat=iostat_val) rho;     if(iostat_val/=0)then;write(*,*)'ERROR: Invalid rho.';stop;end if
    read(*,*,iostat=iostat_val) mu;      if(iostat_val/=0)then;write(*,*)'ERROR: Invalid mu.';stop;end if
    read(*,*,iostat=iostat_val) eps_in;  if(iostat_val/=0)then;write(*,*)'ERROR: Invalid epsilon.';stop;end if
    read(*,*,iostat=iostat_val) L_in;    if(iostat_val/=0)then;write(*,*)'ERROR: Invalid L_scale.';stop;end if

    if(up<=0d0.and.vp<=0d0.and.wp<=0d0)then;write(*,*)'ERROR: At least one fluctuation must be > 0.';stop;end if
    if(U_mean<=0d0)then;write(*,*)'ERROR: U_mean must be > 0.';stop;end if
    if(rho<=0d0)then;write(*,*)'ERROR: rho must be > 0.';stop;end if
    if(mu<=0d0)then;write(*,*)'ERROR: mu must be > 0.';stop;end if
    if(delta<=0d0)then;write(*,*)'ERROR: delta must be > 0.';stop;end if

    nu = mu / rho

    ! ---- TKE ---------------------------------------------------------
    isotropic_flag = .false.
    if (vp <= 0d0 .and. wp <= 0d0) then
        isotropic_flag = .true.
        vp = up;  wp = up
    end if
    k_val = 0.5d0 * (up**2 + vp**2 + wp**2)
    rho_k = rho * k_val
    u_rms = sqrt(2.0d0 * k_val / 3.0d0)
    TI    = u_rms / U_mean
    TI_pct = TI * 100.0d0

    ! ---- Length scale & dissipation ----------------------------------
    if (L_in > 0d0) then
        L_val = L_in
    else
        L_val = 0.09d0 * delta
    end if

    if (eps_in > 0d0) then
        eps_val = eps_in
    else
        eps_val = C_mu**0.75d0 * k_val**1.5d0 / L_val
    end if

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

    ! ---- Taylor microscale -------------------------------------------
    if (eps_val > 0d0) then
        lambda_val = sqrt(10.0d0 * nu * k_val / eps_val)
    else
        lambda_val = 0d0
    end if

    ! ---- Turbulent viscosity -----------------------------------------
    if (eps_val > 0d0) then
        nu_t_ke = C_mu * k_val**2 / eps_val
        omega_val = eps_val / (C_mu * k_val)
        nu_t_kw = k_val / omega_val
    else
        nu_t_ke = 0d0; omega_val = 0d0; nu_t_kw = 0d0
    end if
    nut_ratio = nu_t_ke / nu

    ! ---- Reynolds numbers --------------------------------------------
    if (eps_val > 0d0 .and. nu > 0d0) then
        Re_t = k_val**2 / (nu * eps_val)
    else
        Re_t = 0d0
    end if
    Re_lam = u_rms * lambda_val / nu

    ! ---- Anisotropy --------------------------------------------------
    if (k_val > 0d0) then
        b11 = up**2 / (2.0d0*k_val) - 1.0d0/3.0d0
        b22 = vp**2 / (2.0d0*k_val) - 1.0d0/3.0d0
        b33 = wp**2 / (2.0d0*k_val) - 1.0d0/3.0d0
    else
        b11=0d0; b22=0d0; b33=0d0
    end if
    II_inv  = 2.0d0*(b11**2 + b22**2 + b33**2)
    III_inv = 3.0d0*b11*b22*b33

    ! ---- Production estimate -----------------------------------------
    P_est = nu_t_ke * (U_mean / delta)**2

    ! ---- output ------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   TURBULENT KINETIC ENERGY (k) CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT FLUCTUATIONS --------------------------------------'
    write(*,'(A,F12.6,A)')  '  u_prime (streamwise)    = ', up, ' m/s'
    write(*,'(A,F12.6,A)')  '  v_prime (cross-stream)  = ', vp, ' m/s'
    write(*,'(A,F12.6,A)')  '  w_prime (spanwise)      = ', wp, ' m/s'
    if (isotropic_flag) then
        write(*,'(A)')      '  ** Isotropic assumed (v_prime=w_prime=u_prime) **'
    end if
    write(*,'(A,F12.4,A)')  '  U_mean                  = ', U_mean, ' m/s'
    write(*,'(A,ES14.6,A)') '  Density (rho)           = ', rho, ' kg/m3'
    write(*,'(A,ES14.6,A)') '  Dynamic Viscosity (mu)  = ', mu, ' Pa.s'
    write(*,'(A,ES14.6,A)') '  Kinematic Visc (nu)     = ', nu, ' m2/s'
    write(*,'(A,ES14.6,A)') '  BL Thickness (delta)    = ', delta, ' m'
    write(*,*)

    write(*,'(A)') '--- TURBULENT KINETIC ENERGY --------------------------------'
    write(*,'(A,ES14.6,A)') '  k = 0.5*(u2+v2+w2)     = ', k_val, ' m2/s2'
    write(*,'(A,ES14.6,A)') '  rho*k (per volume)      = ', rho_k, ' Pa'
    write(*,'(A,F12.6,A)')  '  u_rms = sqrt(2k/3)      = ', u_rms, ' m/s'
    write(*,*)

    write(*,'(A)') '--- TURBULENCE INTENSITY ------------------------------------'
    write(*,'(A,F12.6)')    '  TI = u_rms / U_mean     = ', TI
    write(*,'(A,F12.4,A)')  '  TI (percent)            = ', TI_pct, ' %'
    write(*,*)

    write(*,'(A)') '--- DISSIPATION AND SCALES ----------------------------------'
    write(*,'(A,ES14.6,A)') '  epsilon (dissipation)   = ', eps_val, ' m2/s3'
    write(*,'(A,ES14.6,A)') '  Integral Scale L        = ', L_val, ' m'
    write(*,'(A,ES14.6,A)') '  Taylor Microscale lam   = ', lambda_val, ' m'
    write(*,'(A,ES14.6,A)') '  Kolmogorov eta          = ', eta, ' m'
    write(*,'(A,ES14.6,A)') '  Kolmogorov tau_eta      = ', tau_eta, ' s'
    write(*,'(A,ES14.6,A)') '  Kolmogorov u_eta        = ', u_eta, ' m/s'
    write(*,*)

    write(*,'(A)') '--- TURBULENT VISCOSITY -------------------------------------'
    write(*,'(A,ES14.6,A)') '  nu_t (k-epsilon model)  = ', nu_t_ke, ' m2/s'
    write(*,'(A,F14.4)')    '  nu_t / nu               = ', nut_ratio
    write(*,'(A,ES14.6,A)') '  omega (specific dissip)  = ', omega_val, ' 1/s'
    write(*,'(A,ES14.6,A)') '  nu_t (k-omega model)    = ', nu_t_kw, ' m2/s'
    write(*,*)

    write(*,'(A)') '--- REYNOLDS NUMBERS ----------------------------------------'
    write(*,'(A,ES14.6)')   '  Re_t = k^2/(nu*eps)     = ', Re_t
    write(*,'(A,ES14.6)')   '  Re_lambda               = ', Re_lam
    write(*,*)

    write(*,'(A)') '--- ANISOTROPY (Lumley) -------------------------------------'
    write(*,'(A,F12.6)')    '  b11 (streamwise)        = ', b11
    write(*,'(A,F12.6)')    '  b22 (cross-stream)      = ', b22
    write(*,'(A,F12.6)')    '  b33 (spanwise)          = ', b33
    write(*,'(A,ES14.6)')   '  II  invariant           = ', II_inv
    write(*,'(A,ES14.6)')   '  III invariant           = ', III_inv
    if (abs(b11) < 0.01d0 .and. abs(b22) < 0.01d0 .and. abs(b33) < 0.01d0) then
        write(*,'(A)')      '  Status: NEARLY ISOTROPIC'
    else
        write(*,'(A)')      '  Status: ANISOTROPIC'
    end if
    write(*,*)

    write(*,'(A)') '--- PRODUCTION ESTIMATE -------------------------------------'
    write(*,'(A,ES14.6,A)') '  P ~ nu_t*(U/delta)^2    = ', P_est, ' m2/s3'
    write(*,'(A,F12.4)')    '  P / epsilon              = ', P_est / eps_val
    write(*,*)

    ! ---- profile sweep -----------------------------------------------
    write(*,'(A)') '--- PROFILE vs TURBULENCE INTENSITY -------------------------'
    write(*,'(A)') '  TI%       k           epsilon      nu_t/nu      Re_t        eta'
    write(*,'(A)') '  ------------------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        ti_cur = dble(i) * 0.75d0        ! TI% from 0.75 to 30
        urms_c = ti_cur / 100.0d0 * U_mean
        k_c    = 1.5d0 * urms_c**2
        eps_c  = C_mu**0.75d0 * k_c**1.5d0 / L_val
        nut_c  = C_mu * k_c**2 / eps_c / nu
        ret_c  = k_c**2 / (nu * eps_c)
        eta_c  = (nu**3 / eps_c)**0.25d0

        write(*,'(F8.2,2X,ES12.4,2X,ES12.4,2X,F12.2,2X,ES12.4,2X,ES12.4)') &
            ti_cur, k_c, eps_c, nut_c, ret_c, eta_c
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  k = 0.5*(u''^2 + v''^2 + w''^2)'
    write(*,'(A)') '  TI = sqrt(2k/3) / U_mean'
    write(*,'(A)') '  nu_t = C_mu * k^2 / epsilon    (k-epsilon)'
    write(*,'(A)') '  eta = (nu^3/epsilon)^(1/4)      (Kolmogorov)'
    write(*,'(A)') '  lambda = sqrt(10*nu*k/epsilon)  (Taylor)'
    write(*,'(A)') '  Re_t = k^2 / (nu*epsilon)'
    write(*,'(A)') '  b_ij = <u_i u_j>/(2k) - d_ij/3 (anisotropy)'
    write(*,'(A)') '============================================================'

end program turbulent_ke
