program blasius_boundary_layer
    implicit none
    double precision :: U_inf, rho, mu, x_pos, T_wall, T_inf
    double precision :: Pr_in, cp_val, k_cond
    double precision :: nu, Re_x, Re_crit, x_trans, Pr_mol
    ! laminar
    double precision :: d_lam, ds_lam, th_lam, shape_H_lam
    double precision :: Cf_lam, Cfa_lam, tau_lam, ut_lam, FD_lam
    double precision :: dT_lam, Nu_lam, Nua_lam, h_lam, q_lam, St_lam
    ! turbulent
    double precision :: d_turb, ds_turb, th_turb, shape_H_turb
    double precision :: Cf_turb, Cfa_turb, tau_turb, ut_turb, FD_turb
    double precision :: Nu_turb, Nua_turb, h_turb, q_turb, St_turb
    character(len=20) :: regime
    integer :: i, n_pts, ios
    double precision :: xc, rc, dl, dt2, dsl, dst, thl, tht, cfl, cft, nul2, nut2

    read(*,*,iostat=ios) U_inf;  if(ios/=0)then;write(*,*)'ERROR: Invalid U_inf.';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) x_pos;  if(ios/=0)then;write(*,*)'ERROR: Invalid x.';stop;end if
    read(*,*,iostat=ios) T_wall; if(ios/=0)then;write(*,*)'ERROR: Invalid T_wall.';stop;end if
    read(*,*,iostat=ios) T_inf;  if(ios/=0)then;write(*,*)'ERROR: Invalid T_inf.';stop;end if
    read(*,*,iostat=ios) Pr_in;  if(ios/=0)then;write(*,*)'ERROR: Invalid Pr.';stop;end if
    read(*,*,iostat=ios) cp_val; if(ios/=0)then;write(*,*)'ERROR: Invalid cp.';stop;end if
    read(*,*,iostat=ios) k_cond; if(ios/=0)then;write(*,*)'ERROR: Invalid k_cond.';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(x_pos<=0d0)then;write(*,*)'ERROR: x must be > 0.';stop;end if
    if(cp_val<=0d0)then;write(*,*)'ERROR: cp must be > 0.';stop;end if
    if(k_cond<=0d0)then;write(*,*)'ERROR: k_cond must be > 0.';stop;end if

    nu = mu / rho
    Re_x = U_inf * x_pos / nu
    Re_crit = 5d5
    x_trans = Re_crit * nu / U_inf

    if (Pr_in > 0d0) then
        Pr_mol = Pr_in
    else
        Pr_mol = mu * cp_val / k_cond
    end if

    if (Re_x < 3d5) then
        regime = 'LAMINAR'
    else if (Re_x > 1d6) then
        regime = 'TURBULENT'
    else
        regime = 'TRANSITIONAL'
    end if

    ! ---- LAMINAR (Blasius) -------------------------------------------
    d_lam  = 5d0 * x_pos / sqrt(Re_x)
    ds_lam = 1.7208d0 * x_pos / sqrt(Re_x)
    th_lam = 0.664d0 * x_pos / sqrt(Re_x)
    shape_H_lam  = ds_lam / th_lam
    Cf_lam = 0.664d0 / sqrt(Re_x)
    Cfa_lam = 1.328d0 / sqrt(Re_x)
    tau_lam = 0.5d0 * Cf_lam * rho * U_inf**2
    ut_lam  = sqrt(tau_lam / rho)
    FD_lam  = 0.5d0 * Cfa_lam * rho * U_inf**2 * x_pos
    dT_lam  = d_lam / Pr_mol**(1d0/3d0)
    Nu_lam  = 0.332d0 * sqrt(Re_x) * Pr_mol**(1d0/3d0)
    Nua_lam = 0.664d0 * sqrt(Re_x) * Pr_mol**(1d0/3d0)
    h_lam   = Nu_lam * k_cond / x_pos
    q_lam   = h_lam * (T_wall - T_inf)
    if (rho*cp_val*U_inf > 0d0) then
        St_lam = h_lam / (rho * cp_val * U_inf)
    else
        St_lam = 0d0
    end if

    ! ---- TURBULENT (1/7 power law) -----------------------------------
    d_turb  = 0.37d0 * x_pos / Re_x**0.2d0
    ds_turb = 0.0463d0 * x_pos / Re_x**0.2d0
    th_turb = 0.0360d0 * x_pos / Re_x**0.2d0
    shape_H_turb  = ds_turb / th_turb
    Cf_turb = 0.0592d0 / Re_x**0.2d0
    Cfa_turb = 0.074d0 / Re_x**0.2d0
    tau_turb = 0.5d0 * Cf_turb * rho * U_inf**2
    ut_turb  = sqrt(tau_turb / rho)
    FD_turb  = 0.5d0 * Cfa_turb * rho * U_inf**2 * x_pos
    Nu_turb  = 0.0296d0 * Re_x**0.8d0 * Pr_mol**(1d0/3d0)
    Nua_turb = 0.037d0 * Re_x**0.8d0 * Pr_mol**(1d0/3d0)
    h_turb   = Nu_turb * k_cond / x_pos
    q_turb   = h_turb * (T_wall - T_inf)
    if (rho*cp_val*U_inf > 0d0) then
        St_turb = h_turb / (rho * cp_val * U_inf)
    else
        St_turb = 0d0
    end if

    ! ---- Output ------------------------------------------------------
    write(*,'(A)') '============================================================'
    write(*,'(A)') '   BLASIUS BOUNDARY LAYER CALCULATOR'
    write(*,'(A)') '============================================================'
    write(*,*)

    write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
    write(*,'(A,F12.4,A)')  '  U_inf                   = ', U_inf, ' m/s'
    write(*,'(A,F12.4,A)')  '  Density (rho)           = ', rho, ' kg/m3'
    write(*,'(A,ES14.6,A)') '  Dyn. Viscosity (mu)     = ', mu, ' Pa.s'
    write(*,'(A,ES14.6,A)') '  Kinematic Visc (nu)     = ', nu, ' m2/s'
    write(*,'(A,F12.6,A)')  '  x (from leading edge)   = ', x_pos, ' m'
    write(*,'(A,F12.2,A)')  '  T_wall                  = ', T_wall, ' K'
    write(*,'(A,F12.2,A)')  '  T_inf                   = ', T_inf, ' K'
    write(*,'(A,F12.4)')    '  Prandtl Number (Pr)     = ', Pr_mol
    write(*,'(A,F12.2,A)')  '  cp                      = ', cp_val, ' J/kg.K'
    write(*,'(A,F12.6,A)')  '  k (conductivity)        = ', k_cond, ' W/m.K'
    write(*,*)

    write(*,'(A)') '--- REYNOLDS NUMBER & REGIME --------------------------------'
    write(*,'(A,ES14.6)')   '  Re_x                    = ', Re_x
    write(*,'(A,ES14.6)')   '  Re_crit (assumed)       = ', Re_crit
    write(*,'(A,F12.6,A)')  '  x_transition            = ', x_trans, ' m'
    write(*,'(A,A)')        '  Flow Regime             = ', trim(regime)
    write(*,*)

    write(*,'(A)') '--- BOUNDARY LAYER (Laminar - Blasius) ----------------------'
    write(*,'(A,F12.6,A)')  '  delta (BL thick)        = ', d_lam*1d3, ' mm'
    write(*,'(A,F12.6,A)')  '  delta* (displacement)   = ', ds_lam*1d3, ' mm'
    write(*,'(A,F12.6,A)')  '  theta (momentum)        = ', th_lam*1d3, ' mm'
    write(*,'(A,F12.6)')    '  H = delta*/theta        = ', shape_H_lam
    write(*,'(A,F12.6,A)')  '  delta_T (thermal)       = ', dT_lam*1d3, ' mm'
    write(*,*)

    write(*,'(A)') '--- BOUNDARY LAYER (Turbulent - 1/7 law) -------------------'
    write(*,'(A,F12.6,A)')  '  delta (BL thick)        = ', d_turb*1d3, ' mm'
    write(*,'(A,F12.6,A)')  '  delta* (displacement)   = ', ds_turb*1d3, ' mm'
    write(*,'(A,F12.6,A)')  '  theta (momentum)        = ', th_turb*1d3, ' mm'
    write(*,'(A,F12.6)')    '  H = delta*/theta        = ', shape_H_turb
    write(*,*)

    write(*,'(A)') '--- SKIN FRICTION -------------------------------------------'
    write(*,'(A,20X,A14,A14)')    '  ','Laminar','Turbulent'
    write(*,'(A,ES14.6,2X,ES14.6)')   '  Cf (local)  ', Cf_lam, Cf_turb
    write(*,'(A,ES14.6,2X,ES14.6)')   '  Cf (avg)    ', Cfa_lam, Cfa_turb
    write(*,'(A,ES14.6,2X,ES14.6,A)') '  tau_w       ', tau_lam, tau_turb, ' Pa'
    write(*,'(A,ES14.6,2X,ES14.6,A)') '  u_tau       ', ut_lam, ut_turb, ' m/s'
    write(*,'(A,ES14.6,2X,ES14.6,A)') '  F_D/width   ', FD_lam, FD_turb, ' N/m'
    write(*,*)

    write(*,'(A)') '--- HEAT TRANSFER -------------------------------------------'
    write(*,'(A,20X,A14,A14)')    '  ','Laminar','Turbulent'
    write(*,'(A,ES14.6,2X,ES14.6)')   '  Nu_x (local)', Nu_lam, Nu_turb
    write(*,'(A,ES14.6,2X,ES14.6)')   '  Nu_x (avg)  ', Nua_lam, Nua_turb
    write(*,'(A,ES14.6,2X,ES14.6,A)') '  h (local)   ', h_lam, h_turb, ' W/m2.K'
    write(*,'(A,ES14.6,2X,ES14.6,A)') '  q (flux)    ', q_lam, q_turb, ' W/m2'
    write(*,'(A,ES14.6,2X,ES14.6)')   '  St          ', St_lam, St_turb
    write(*,*)

    ! ---- Profile sweep -----------------------------------------------
    write(*,'(A)') '--- PROFILE vs x -------------------------------------------'
    write(*,'(A)') '  x[m]      Re_x        d_lam[mm]   d_turb[mm]  Cf_lam      Cf_turb     Nu_lam'
    write(*,'(A)') '  ----------------------------------------------------------------------------------'

    n_pts = 40
    do i = 1, n_pts
        xc = 0.01d0 + (max(x_pos*2d0, 1d0) - 0.01d0) * dble(i-1) / dble(n_pts-1)
        rc = U_inf * xc / nu
        if (rc < 1d0) rc = 1d0
        dl  = 5d0 * xc / sqrt(rc) * 1d3
        dt2 = 0.37d0 * xc / rc**0.2d0 * 1d3
        cfl = 0.664d0 / sqrt(rc)
        cft = 0.0592d0 / rc**0.2d0
        nul2 = 0.332d0 * sqrt(rc) * Pr_mol**(1d0/3d0)

        write(*,'(F8.4,2X,ES12.4,2X,F10.4,2X,F10.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            xc, rc, dl, dt2, cfl, cft, nul2
    end do

    write(*,*)
    write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
    write(*,'(A)') '  Blasius:  delta = 5x/sqrt(Re_x)'
    write(*,'(A)') '  Blasius:  Cf = 0.664/sqrt(Re_x)'
    write(*,'(A)') '  Turb:    delta = 0.37x/Re_x^(1/5)'
    write(*,'(A)') '  Turb:    Cf = 0.0592/Re_x^(1/5)'
    write(*,'(A)') '  Nu_lam = 0.332*Re_x^0.5*Pr^(1/3)'
    write(*,'(A)') '  Nu_turb = 0.0296*Re_x^0.8*Pr^(1/3)'
    write(*,'(A)') '============================================================'

end program blasius_boundary_layer
