Reflux Partial Condenser Sizer (Dephlegmator)

program reflux_partial_condenser_dephlegmator
    implicit none
    integer :: iostat_val, mixture_type, num_tubes
    double precision :: m_dot_vapor_kgh, Tv_in_C, yin_mole_pct, T_cool_in_C, reflux_fraction_target
    double precision :: D_tube_mm, L_tube_m, m_dot_cool_kgh
    double precision :: rho_L, rho_v, cp_v, mu_v, k_v, hfg_kJkg, sigma_Nm
    double precision :: m_v_kgs, m_cond_kgs, m_over_kgs, D_i_m, A_flow_tot, V_vapor_ms, V_flood_ms
    double precision :: h_cond, h_gas, h_cool, Z_ratio, U_overall, A_total_m2, Q_duty_kW
    double precision :: y_overhead_mole_pct, flooding_margin_pct
    character(len=32) :: flooding_status
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: G_ACC = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) mixture_type           ! 1=Ethanol-Water, 2=Methanol-Water, 3=Benzene-Toluene, 4=Acetone-Water
    read(*,*,iostat=iostat_val) m_dot_vapor_kgh        ! Vapor Feed Flow from Column [kg/h] (e.g. 1500.0)
    read(*,*,iostat=iostat_val) Tv_in_C                ! Vapor Feed Temperature [deg C] (e.g. 84.0)
    read(*,*,iostat=iostat_val) yin_mole_pct           ! Light Component Inlet Mole Fraction [%] (e.g. 65.0)
    read(*,*,iostat=iostat_val) reflux_fraction_target ! Condensed Reflux Split Fraction [0 to 0.95] (e.g. 0.60)
    read(*,*,iostat=iostat_val) T_cool_in_C            ! Shell Cooling Water Inlet Temp [deg C] (e.g. 25.0)
    read(*,*,iostat=iostat_val) num_tubes              ! Number of Vertical Condenser Tubes (e.g. 80)
    read(*,*,iostat=iostat_val) D_tube_mm              ! Tube Inner Diameter [mm] (e.g. 32.0)
    read(*,*,iostat=iostat_val) L_tube_m               ! Tube Active Length [m] (e.g. 3.0)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for reflux condenser calculation.'
        stop
    end if

    if (m_dot_vapor_kgh <= 0.0d0 .or. num_tubes <= 0 .or. D_tube_mm <= 0.0d0 .or. L_tube_m <= 0.0d0) then
        write(*,*) 'ERROR: Flow rate, tubes count, and tube dimensions must be positive.'
        stop
    end if

    ! Mixture properties
    if (mixture_type == 1) then ! Ethanol-Water
        rho_L = 820.0d0; rho_v = 1.65d0; cp_v = 1900.0d0; mu_v = 0.000012d0; k_v = 0.024d0; hfg_kJkg = 1200.0d0; sigma_Nm = 0.022d0
    else if (mixture_type == 2) then ! Methanol-Water
        rho_L = 840.0d0; rho_v = 1.35d0; cp_v = 1850.0d0; mu_v = 0.000011d0; k_v = 0.022d0; hfg_kJkg = 1350.0d0; sigma_Nm = 0.024d0
    else if (mixture_type == 3) then ! Benzene-Toluene
        rho_L = 830.0d0; rho_v = 2.80d0; cp_v = 1450.0d0; mu_v = 0.000010d0; k_v = 0.018d0; hfg_kJkg = 390.0d0;  sigma_Nm = 0.025d0
    else ! Acetone-Water
        rho_L = 850.0d0; rho_v = 2.10d0; cp_v = 1600.0d0; mu_v = 0.000011d0; k_v = 0.020d0; hfg_kJkg = 850.0d0;  sigma_Nm = 0.023d0
    end if

    m_v_kgs = m_dot_vapor_kgh / 3600.0d0
    m_cond_kgs = m_v_kgs * reflux_fraction_target
    m_over_kgs = m_v_kgs - m_cond_kgs

    D_i_m = D_tube_mm / 1000.0d0
    A_flow_tot = dble(num_tubes) * PI * (D_i_m**2) / 4.0d0
    A_total_m2 = dble(num_tubes) * PI * D_i_m * L_tube_m

    V_vapor_ms = m_v_kgs / (rho_v * max(1.0d-6, A_flow_tot))

    ! Hewitt & Wallis Flooding Limit Velocity [m/s]
    V_flood_ms = 0.85d0 * sqrt(G_ACC * D_i_m * (rho_L - rho_v) / rho_v) * &
                 ((1.0d0 + ((rho_v/rho_L)**0.25d0) * sqrt(max(0.01d0, m_cond_kgs/m_v_kgs)))**(-2.0d0))

    flooding_margin_pct = (V_vapor_ms / max(0.1d0, V_flood_ms)) * 100.0d0

    if (flooding_margin_pct > 100.0d0) then
        flooding_status = 'WARNING: COLUMN FLOODED / BLOWOUT'
    else if (flooding_margin_pct > 80.0d0) then
        flooding_status = 'HIGH LIQUID HOLDUP (NEAR FLOODING)'
    else
        flooding_status = 'SAFE COUNTER-CURRENT FLOW'
    end if

    ! Thermal sizing (Silver-Bell-Ghaly model)
    h_cond = 2800.0d0
    h_gas = 85.0d0
    h_cool = 4500.0d0
    Z_ratio = 0.15d0 ! Sensible gas cooling fraction approx 15%

    U_overall = 1.0d0 / ((1.0d0/h_cond) + (Z_ratio/h_gas) + (0.002d0/16.0d0) + (1.0d0/h_cool))

    Q_duty_kW = (m_cond_kgs * hfg_kJkg) + (m_v_kgs * cp_v * 8.0d0 / 1000.0d0)

    ! Vapor enrichment (Rayleigh dephlegmation fractionation boost)
    y_overhead_mole_pct = min(98.5d0, yin_mole_pct + 12.0d0 * reflux_fraction_target)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — REFLUX CONDENSER (DEPHLEGMATOR SIZING)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Vapor Velocity / Flood Limit = ', V_vapor_ms, ' m/s / ', V_flood_ms, ' m/s'
    write(*,'(A,F10.1,A,A)')       'Flooding Approach Fraction   = ', flooding_margin_pct, ' % | ', trim(flooding_status)
    write(*,'(A,F10.1,A,F10.1,A)') 'Inlet Mole / Overhead Mole   = ', yin_mole_pct, ' % / ', y_overhead_mole_pct, ' % (Enriched)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'CONDENSER THERMAL DUTY (Q)   = ', Q_duty_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Overall Heat Transfer Coeff U= ', U_overall, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'Reflux Condensate Return     = ', m_cond_kgs*3600.0d0, ' kg/h'
    write(*,'(A,F10.1,A)')  'Purified Overhead Vapor Flow = ', m_over_kgs*3600.0d0, ' kg/h'
    write(*,'(A,F10.2,A,I4,A)') 'Total Surface Area       = ', A_total_m2, ' m2 (', num_tubes, ' vertical tubes)'
    write(*,'(A)') '============================================================'

end program reflux_partial_condenser_dephlegmator