program entry_length
    implicit none
    double precision :: D, V, Ti, Tw, Lp, rho, mu, kf, Pr, cp, bc_val
    integer :: fl_type, bc_type
    
    double precision :: Tf, dT, nu, Re, Gz, Lh, Lt, Nufd, Nud, h, As, Q, To, g, f, pi
    integer :: iostat_val, i
    double precision :: temp_val(6)
    double precision :: xs, Gzs, Nuds, hs
    
    pi = 3.141592653589793d0
    
    ! Read inputs sequentially
    read(*,*,iostat=iostat_val) D
    read(*,*,iostat=iostat_val) V
    read(*,*,iostat=iostat_val) Ti
    read(*,*,iostat=iostat_val) Tw
    read(*,*,iostat=iostat_val) Lp
    read(*,*,iostat=iostat_val) fl_type
    
    ! Read custom parameters
    read(*,*,iostat=iostat_val) temp_val(1) ! rho
    read(*,*,iostat=iostat_val) temp_val(2) ! mu
    read(*,*,iostat=iostat_val) temp_val(3) ! kf
    read(*,*,iostat=iostat_val) temp_val(4) ! Pr
    read(*,*,iostat=iostat_val) temp_val(5) ! cp
    read(*,*,iostat=iostat_val) bc_val      ! boundary condition 1=UWT, 2=UHF
    
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all input values.'
        stop
    end if
    
    bc_type = int(bc_val)
    Tf = (Ti + Tw) / 2.0d0
    dT = abs(Tw - Ti)
    
    ! Select defaults based on fluid type
    if (fl_type == 1) then
        ! Air
        rho = 1.177d0
        mu = 1.85d-5
        kf = 0.0263d0
        Pr = 0.71d0
        cp = 1007d0
    else if (fl_type == 2) then
        ! Water
        rho = 997d0
        mu = 8.9d-4
        kf = 0.613d0
        Pr = 6.13d0
        cp = 4180d0
    else
        ! Engine Oil
        rho = 870d0
        mu = 0.05d0
        kf = 0.14d0
        Pr = 500d0
        cp = 2000d0
    end if
    
    ! Override with custom parameters if provided
    if (temp_val(1) > 0.0d0) rho = temp_val(1)
    if (temp_val(2) > 0.0d0) mu = temp_val(2)
    if (temp_val(3) > 0.0d0) kf = temp_val(3)
    if (temp_val(4) > 0.0d0) Pr = temp_val(4)
    if (temp_val(5) > 0.0d0) cp = temp_val(5)
    
    nu = mu / rho
    Re = rho * V * D / mu
    Gz = (D / Lp) * Re * Pr
    
    ! Hydrodynamic and thermal entry lengths
    if (Re <= 2300.0d0) then
        ! Laminar
        Lh = 0.05d0 * Re * D
        Lt = Lh * Pr
    else
        ! Turbulent
        Lh = 10.0d0 * D
        Lt = 10.0d0 * D
    end if
    
    ! Nusselt calculations
    if (Re <= 2300.0d0) then
        ! Laminar flow
        if (bc_type == 1) then
            ! Constant Wall Temp (UWT)
            Nufd = 3.66d0
            Nud = 3.66d0 + 0.0668d0 * Gz / (1.0d0 + 0.04d0 * Gz**(2.0d0/3.0d0))
        else
            ! Constant Heat Flux (UHF)
            Nufd = 4.36d0
            if (Gz > 0.001d0) then
                Nud = 1.953d0 * Gz**(1.0d0/3.0d0)
                if (Nud < 4.36d0) Nud = 4.36d0
            else
                Nud = 4.36d0
            end if
        end if
    else
        ! Turbulent flow (Gnielinski + Entry length correction)
        f = (0.790d0 * log(Re) - 1.64d0)**(-2)
        Nufd = (f / 8.0d0) * (Re - 1000.0d0) * Pr / &
               (1.0d0 + 12.7d0 * sqrt(f / 8.0d0) * (Pr**(2.0d0/3.0d0) - 1.0d0))
        Nud = Nufd * (1.0d0 + (D / Lp)**0.7d0)
    end if
    
    h = Nud * kf / D
    As = pi * D * Lp
    
    ! Heat transfer and outlet temperature
    if (bc_type == 1) then
        ! UWT
        To = Tw - (Tw - Ti) * exp(-h * As / (rho * V * (pi * D**2 / 4.0d0) * cp))
        Q = rho * V * (pi * D**2 / 4.0d0) * cp * (To - Ti)
    else
        ! UHF
        Q = h * As * dT
        To = Ti + Q / (rho * V * (pi * D**2 / 4.0d0) * cp)
    end if
    
    write(*,'(A)') '============================================'
    write(*,'(A)') '  INTERNAL FLOW ENTRY LENGTH HEAT TRANSFER'
    write(*,'(A)') '============================================'
    write(*,*)
    write(*,'(A,ES14.4)') '  Reynolds Re             = ', Re
    write(*,'(A,F12.4)')    '  Graetz Number           = ', Gz
    write(*,'(A,F12.4)')    '  Nu Developing           = ', Nud
    write(*,'(A,F12.4)')    '  Nu Fully                = ', Nufd
    write(*,'(A,F12.4,A)')  '  Coeff h                 = ', h, ' W/m2K'
    write(*,'(A,F12.4,A)')  '  Transfer Q              = ', abs(Q), ' W'
    write(*,'(A,F12.2,A)')  '  T_out                   = ', To, ' C'
    write(*,'(A,F12.4,A)')  '  Entry Len L_h           = ', Lh, ' m'
    write(*,'(A,F12.4,A)')  '  Entry Len L_t           = ', Lt, ' m'
    write(*,*)
    
    write(*,'(A)') '--- LENGTH SWEEP ---'
    write(*,'(A)') '  x[m]       Gz           Nu_local     h_local[W/m2K]'
    write(*,'(A)') '  --------------------------------------------------------'
    do i=1,25
        xs = max(0.001d0 * Lp, 1.0d-4) + (Lp - max(0.001d0 * Lp, 1.0d-4))*dble(i-1)/24.0d0
        Gzs = (D / xs) * Re * Pr
        if (Re <= 2300.0d0) then
            if (bc_type == 1) then
                Nuds = 3.66d0 + 0.0668d0 * Gzs / (1.0d0 + 0.04d0 * Gzs**(2.0d0/3.0d0))
            else
                if (Gzs > 0.001d0) then
                    Nuds = 1.953d0 * Gzs**(1.0d0/3.0d0)
                    if (Nuds < 4.36d0) Nuds = 4.36d0
                else
                    Nuds = 4.36d0
                end if
            end if
        else
            Nuds = Nufd * (1.0d0 + (D / xs)**0.7d0)
        end if
        hs = Nuds * kf / D
        write(*,'(2X,F8.4,2X,ES10.3,2X,F12.4,2X,F14.4)') xs, Gzs, Nuds, hs
    end do
    write(*,*)
end program entry_length