program two_phase_lockhart_martinelli
    implicit none
    integer :: iostat_val, pipe_orient
    double precision :: m_dot_tot_kgs, vapor_qual_x, D_pipe_mm, L_pipe_m, P_sat_bar
    double precision :: rho_L, rho_G, mu_L, mu_G, g, D_m, A_pipe, m_dot_L, m_dot_G
    double precision :: G_mass_flux, Re_L, Re_G, f_L, f_G, dp_L_Pam, dp_G_Pam
    double precision :: X_LM, C_chisholm, phi2_L, dp_tp_Pam, dp_tot_kPa, void_frac, liquid_holdup
    character(len=32) :: flow_regime
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) pipe_orient     ! 1=Horizontal, 2=Vertical Upward
    read(*,*,iostat=iostat_val) m_dot_tot_kgs   ! Total Mass Flow Rate [kg/s] (e.g. 5.0)
    read(*,*,iostat=iostat_val) vapor_qual_x    ! Vapor Quality x [0 to 1] (e.g. 0.25)
    read(*,*,iostat=iostat_val) D_pipe_mm       ! Pipe Inner Diameter [mm] (e.g. 80.0)
    read(*,*,iostat=iostat_val) L_pipe_m        ! Pipe Length [m] (e.g. 50.0)
    read(*,*,iostat=iostat_val) rho_L           ! Liquid Density [kg/m3] (e.g. 998.0)
    read(*,*,iostat=iostat_val) rho_G           ! Gas/Vapor Density [kg/m3] (e.g. 2.4)
    read(*,*,iostat=iostat_val) mu_L            ! Liquid Viscosity [Pa.s] (e.g. 0.0010)
    read(*,*,iostat=iostat_val) mu_G            ! Gas Viscosity [Pa.s] (e.g. 0.000018)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for two-phase Lockhart-Martinelli calculation.'
        stop
    end if

    if (D_pipe_mm <= 0.0d0 .or. m_dot_tot_kgs <= 0.0d0 .or. vapor_qual_x < 0.0d0 .or. vapor_qual_x > 1.0d0) then
        write(*,*) 'ERROR: Pipe diameter, flow must be positive and quality x in [0, 1].'
        stop
    end if

    g = 9.80665d0
    D_m = D_pipe_mm / 1000.0d0
    A_pipe = PI * (D_m**2) / 4.0d0
    G_mass_flux = m_dot_tot_kgs / A_pipe

    m_dot_L = m_dot_tot_kgs * (1.0d0 - vapor_qual_x)
    m_dot_G = m_dot_tot_kgs * vapor_qual_x

    ! Superficial Reynolds numbers
    Re_L = (4.0d0 * m_dot_L) / (PI * D_m * max(1.0d-6, mu_L))
    Re_G = (4.0d0 * m_dot_G) / (PI * D_m * max(1.0d-6, mu_G))

    ! Liquid friction
    if (Re_L < 2000.0d0) then
        f_L = 64.0d0 / max(1.0d0, Re_L)
    else
        f_L = 0.3164d0 / (Re_L**0.25d0)
    end if
    dp_L_Pam = f_L * (1.0d0 / D_m) * ((m_dot_L / A_pipe)**2) / (2.0d0 * rho_L)

    ! Gas friction
    if (Re_G < 2000.0d0) then
        f_G = 64.0d0 / max(1.0d0, Re_G)
    else
        f_G = 0.3164d0 / (Re_G**0.25d0)
    end if
    dp_G_Pam = f_G * (1.0d0 / D_m) * ((m_dot_G / A_pipe)**2) / (2.0d0 * rho_G)

    ! Lockhart-Martinelli parameter X
    if (dp_G_Pam > 1.0d-9) then
        X_LM = sqrt(dp_L_Pam / dp_G_Pam)
    else
        X_LM = 1000.0d0
    end if

    ! Chisholm C factor
    if (Re_L >= 2000.0d0 .and. Re_G >= 2000.0d0) then
        C_chisholm = 20.0d0 ! tt (turbulent liquid, turbulent gas)
    else if (Re_L < 2000.0d0 .and. Re_G >= 2000.0d0) then
        C_chisholm = 12.0d0 ! vt (viscous liquid, turbulent gas)
    else if (Re_L >= 2000.0d0 .and. Re_G < 2000.0d0) then
        C_chisholm = 10.0d0 ! tv (turbulent liquid, viscous gas)
    else
        C_chisholm = 5.0d0  ! vv (viscous liquid, viscous gas)
    end if

    phi2_L = 1.0d0 + (C_chisholm / max(0.001d0, X_LM)) + (1.0d0 / max(0.0001d0, X_LM**2))
    dp_tp_Pam = phi2_L * dp_L_Pam
    dp_tot_kPa = (dp_tp_Pam * L_pipe_m) / 1000.0d0

    ! Liquid holdup & void fraction (Butterworth correlation)
    liquid_holdup = 1.0d0 / sqrt(1.0d0 + (21.0d0 / max(0.01d0, X_LM)) + (1.0d0 / max(0.0001d0, X_LM**2)))
    void_frac = 1.0d0 - liquid_holdup

    ! Flow regime classification (Baker & Taitel-Dukler)
    if (vapor_qual_x < 0.05d0) then
        flow_regime = 'BUBBLY / DISPERSED'
    else if (vapor_qual_x < 0.35d0 .and. X_LM > 2.0d0) then
        flow_regime = 'SLUG / INTERMITTENT'
    else if (vapor_qual_x < 0.35d0 .and. X_LM <= 2.0d0) then
        flow_regime = 'STRATIFIED-WAVY'
    else if (vapor_qual_x < 0.85d0) then
        flow_regime = 'ANNULAR (LIQUID FILM)'
    else
        flow_regime = 'MIST / DROPLET ENTRAINED'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TWO-PHASE FLOW & LOCKHART-MARTINELLI'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.3)')   'Mass Flow / Vapor Quality x  = ', m_dot_tot_kgs, ' kg/s / ', vapor_qual_x
    write(*,'(A,F10.1,A,F10.1,A)') 'Pipe ID / Pipe Length        = ', D_pipe_mm, ' mm / ', L_pipe_m, ' m'
    write(*,'(A,F10.3,A,F10.1)')   'Lockhart-Martinelli X / C    = ', X_LM, ' / ', C_chisholm
    write(*,'(A,F10.2)')           'Two-Phase Multiplier (phi2_L)= ', phi2_L
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'TWO-PHASE PRESSURE DROP    = ', dp_tot_kPa, ' kPa'
    write(*,'(A,F10.2,A)')  'Frictional Gradient        = ', dp_tp_Pam, ' Pa/m'
    write(*,'(A,F10.1,A,F10.1,A)') 'Void Fraction / Holdup     = ', void_frac*100.0d0, ' % / ', liquid_holdup*100.0d0, ' %'
    write(*,'(A,A)')        'PREDICTED FLOW REGIME      = ', trim(flow_regime)
    write(*,'(A)') '============================================================'

end program two_phase_lockhart_martinelli