program prandtl_mixing_length
    implicit none
    double precision, parameter :: PI = 3.141592653589793d0

    double precision :: y, delta, U_inf, rho, mu, kappa, A_plus, dp_dx
    double precision :: nu, Re_delta, Cf, tau_w, u_tau, delta_nu
    double precision :: y_plus, u_plus, u_local, du_dy
    double precision :: l_pr, l_vd, l_out, l_comb
    double precision :: nu_t, tau_turb, nut_ratio
    character(len=30) :: layer
    integer :: i, n_pts, iostat_val
    double precision :: eta, yc, ypc, upc, ldc, nrc, bc, duy

    ! ---- read inputs --------------------------------------------------
    read(*,*,iostat=iostat_val) y;       if(iostat_val/=0)then;write(*,*)'ERROR: Invalid y.';stop;end if
    read(*,*,iostat=iostat_val) delta;   if(iostat_val/=0)then;write(*,*)'ERROR: Invalid delta.';stop;end if
    read(*,*,iostat=iostat_val) U_inf;   if(iostat_val/=0)then;write(*,*)'ERROR: Invalid U_inf.';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) kappa;   if(iostat_val/=0)then;write(*,*)'ERROR: Invalid kappa.';stop;end if
    read(*,*,iostat=iostat_val) A_plus;  if(iostat_val/=0)then;write(*,*)'ERROR: Invalid A_plus.';stop;end if
    read(*,*,iostat=iostat_val) dp_dx;   if(iostat_val/=0)then;write(*,*)'ERROR: Invalid dp_dx.';stop;end if

    if(y<=0d0)then;write(*,*)'ERROR: y must be > 0.';stop;end if
    if(delta<=0d0)then;write(*,*)'ERROR: delta must be > 0.';stop;end if
    if(U_inf<=0d0)then;write(*,*)'ERROR: U_inf 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(kappa<=0d0)then;write(*,*)'ERROR: kappa must be > 0.';stop;end if
    if(A_plus<=0d0)then;write(*,*)'ERROR: A_plus must be > 0.';stop;end if

    nu = mu / rho
    Re_delta = rho * U_inf * delta / mu

    ! Skin friction (power-law estimate for turbulent BL)
    if (Re_delta > 1.0d0) then
        Cf = 0.027d0 * Re_delta**(-1.0d0/7.0d0)
    else
        Cf = 0.01d0
    end if
    tau_w   = 0.5d0 * Cf * rho * U_inf**2
    u_tau   = sqrt(tau_w / rho)
    delta_nu = nu / u_tau

    ! Wall-units at y
    y_plus = y * u_tau / nu

    ! Velocity from law of the wall
    if (y_plus <= 5.0d0) then
        u_plus = y_plus
        layer  = 'VISCOUS SUBLAYER'
    else if (y_plus <= 30.0d0) then
        u_plus = 5.0d0 * log(y_plus) - 3.05d0
        layer  = 'BUFFER LAYER'
    else
        u_plus = (1.0d0/kappa)*log(y_plus) + 5.0d0
        layer  = 'LOG-LAW REGION'
    end if
    if (y > delta) then
        u_plus = U_inf / u_tau
        layer  = 'OUTER / WAKE'
    end if
    u_local = u_plus * u_tau

    ! Mixing lengths
    l_pr  = kappa * y                                         ! Prandtl original
    l_vd  = kappa * y * (1.0d0 - exp(-y_plus / A_plus))      ! Van Driest
    l_out = 0.09d0 * delta                                    ! Clauser outer
    l_comb = min(l_vd, l_out)                                 ! combined

    ! du/dy from law of the wall
    if (y_plus <= 5.0d0) then
        du_dy = u_tau**2 / nu           ! du+/dy+ = 1  =>  du/dy = u_tau / delta_nu
    else if (y_plus <= 30.0d0) then
        du_dy = 5.0d0 * u_tau / y       ! du+/dy+ = 5/y+
    else
        du_dy = u_tau / (kappa * y)      ! du+/dy+ = 1/(kappa*y+)
    end if

    ! Turbulent quantities
    nu_t      = l_comb**2 * abs(du_dy)
    tau_turb  = rho * l_comb**2 * du_dy**2
    nut_ratio = nu_t / nu

    ! ---- output -------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   PRANDTL MIXING LENGTH CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- FLOW CONDITIONS -----------------------------------------'
    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,F12.4,A)')  '  Freestream Velocity     = ', U_inf, ' m/s'
    write(*,'(A,ES14.6,A)') '  BL Thickness (delta)    = ', delta, ' m'
    write(*,'(A,ES14.6)')   '  Re_delta                = ', Re_delta
    write(*,'(A,F12.6)')    '  Von Karman kappa        = ', kappa
    write(*,'(A,F12.4)')    '  Van Driest A+           = ', A_plus
    write(*,'(A,ES14.6,A)') '  Pressure Gradient       = ', dp_dx, ' Pa/m'
    write(*,*)

    write(*,'(A)') '--- WALL PARAMETERS -----------------------------------------'
    write(*,'(A,ES14.6)')   '  Skin Friction Cf        = ', Cf
    write(*,'(A,ES14.6,A)') '  Wall Shear Stress tau_w = ', tau_w, ' Pa'
    write(*,'(A,F12.6,A)')  '  Friction Velocity u_tau = ', u_tau, ' m/s'
    write(*,'(A,ES14.6,A)') '  Viscous Length delta_nu  = ', delta_nu, ' m'
    write(*,*)

    write(*,'(A)') '--- WALL UNITS AT y -----------------------------------------'
    write(*,'(A,ES14.6,A)') '  y (position)            = ', y, ' m'
    write(*,'(A,F12.4)')    '  y/delta                 = ', y/delta
    write(*,'(A,F14.4)')    '  y+                      = ', y_plus
    write(*,'(A,F14.6)')    '  u+                      = ', u_plus
    write(*,'(A,F12.4,A)')  '  u (local velocity)      = ', u_local, ' m/s'
    write(*,'(A,A)')        '  Layer Classification    = ', trim(layer)
    write(*,*)

    write(*,'(A)') '--- MIXING LENGTH -------------------------------------------'
    write(*,'(A,ES14.6,A)') '  l (Prandtl kappa*y)     = ', l_pr, ' m'
    write(*,'(A,ES14.6,A)') '  l (Van Driest damped)   = ', l_vd, ' m'
    write(*,'(A,ES14.6,A)') '  l (outer Clauser)       = ', l_out, ' m'
    write(*,'(A,ES14.6,A)') '  l (combined min)        = ', l_comb, ' m'
    write(*,'(A,F12.6)')    '  l_combined / delta      = ', l_comb / delta
    write(*,'(A,F12.6)')    '  l_combined / y          = ', l_comb / y
    write(*,*)

    write(*,'(A)') '--- TURBULENT QUANTITIES ------------------------------------'
    write(*,'(A,ES14.6,A)') '  du/dy (velocity grad)   = ', du_dy, ' 1/s'
    write(*,'(A,ES14.6,A)') '  nu_t (eddy viscosity)   = ', nu_t, ' m2/s'
    write(*,'(A,F14.4)')    '  nu_t / nu               = ', nut_ratio
    write(*,'(A,ES14.6,A)') '  tau_turb (Reynolds str)  = ', tau_turb, ' Pa'
    write(*,'(A,F12.6)')    '  tau_turb / tau_w        = ', tau_turb / tau_w
    write(*,*)

    ! ---- profile sweep ------------------------------------------------
    write(*,'(A)') '--- PROFILE vs y/delta --------------------------------------'
    write(*,'(A)') '  y/delta     y+           u+           l/delta      nu_t/nu'
    write(*,'(A)') '  ----------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        eta = dble(i) / dble(n_pts)            ! y/delta = 0.025 .. 1.0
        yc  = eta * delta
        ypc = yc * u_tau / nu
        ! velocity
        if (ypc <= 5.0d0) then
            upc = ypc
        else if (ypc <= 30.0d0) then
            upc = 5.0d0 * log(ypc) - 3.05d0
        else
            upc = (1.0d0/kappa)*log(ypc) + 5.0d0
        end if
        ! mixing length
        bc  = kappa * yc * (1.0d0 - exp(-ypc / A_plus))
        ldc = min(bc, 0.09d0*delta) / delta
        ! du/dy
        if (ypc <= 5.0d0) then
            duy = u_tau**2 / nu
        else if (ypc <= 30.0d0) then
            duy = 5.0d0 * u_tau / yc
        else
            duy = u_tau / (kappa * yc)
        end if
        nrc = (min(bc, 0.09d0*delta))**2 * abs(duy) / nu

        write(*,'(F8.4,2X,F12.2,2X,F12.4,2X,F12.6,2X,F12.2)') &
            eta, ypc, upc, ldc, nrc
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Prandtl:   l = kappa * y'
    write(*,'(A)') '  Van Driest: l = kappa*y*[1-exp(-y+/A+)]'
    write(*,'(A)') '  Outer:     l = 0.09 * delta  (Clauser)'
    write(*,'(A)') '  Combined:  l = min(l_VD, l_outer)'
    write(*,'(A)') '  tau_turb = rho * l^2 * (du/dy)^2'
    write(*,'(A)') '  nu_t = l^2 * |du/dy|'
    write(*,'(A)') '  Law of wall: u+ = y+ (visc), u+ = (1/k)ln(y+)+5 (log)'
    write(*,'(A)') '============================================================'

end program prandtl_mixing_length