program flow_boiling_critical_heat_flux
    implicit none
    integer :: iostat_val
    double precision :: D_tube_mm, L_heated_m, G_massflux, P_bar, DeltaT_sub_C, q_actual_kWm2
    double precision :: D_m, A_flow_m2, A_heat_m2, Tsat_C, hfg_kJkg, rho_L, rho_V
    double precision :: Delta_h_sub_kJkg, q_CHF_kWm2, DNBR_margin, Q_crit_kW, x_crit_exit
    double precision :: p_MPa, n_pow, A_param, C_param
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) D_tube_mm      ! Inside Diameter [mm] (e.g. 10.0)
    read(*,*,iostat=iostat_val) L_heated_m     ! Heated Length [m] (e.g. 2.0)
    read(*,*,iostat=iostat_val) G_massflux     ! Mass Flux G [kg/(m2.s)] (e.g. 2000.0)
    read(*,*,iostat=iostat_val) P_bar          ! System Pressure [bar] (e.g. 70.0)
    read(*,*,iostat=iostat_val) DeltaT_sub_C   ! Inlet Subcooling [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) q_actual_kWm2  ! Actual Applied Heat Flux [kW/m2] (e.g. 850.0)

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

    if (D_tube_mm <= 0.0d0 .or. L_heated_m <= 0.0d0 .or. G_massflux <= 0.0d0 .or. P_bar <= 0.0d0) then
        write(*,*) 'ERROR: Tube geometry, mass flux, and pressure must be positive.'
        stop
    end if

    D_m = D_tube_mm * 1.0d-3
    A_flow_m2 = (PI / 4.0d0) * (D_m**2)
    A_heat_m2 = PI * D_m * L_heated_m

    p_MPa = P_bar * 0.1d0
    ! Approximate steam saturation properties
    Tsat_C = 100.0d0 + 28.0d0 * (P_bar**0.25d0) * 1.8d0 ! approx saturation curve
    if (P_bar > 1.0d0) Tsat_C = 42.6776d0 + 179.916d0 * (P_bar**0.18d0) - 273.15d0 + 150.0d0
    if (P_bar >= 70.0d0) Tsat_C = 285.8d0
    if (P_bar <= 1.013d0) Tsat_C = 100.0d0

    hfg_kJkg = max(800.0d0, 2257.0d0 - 14.5d0 * P_bar)
    rho_L = max(600.0d0, 1000.0d0 - 2.5d0 * P_bar)
    rho_V = 0.59d0 * (P_bar**0.95d0)

    Delta_h_sub_kJkg = 4.2d0 * DeltaT_sub_C ! cp_L * DeltaT_sub

    ! Bowring / Katto-Ohno Flow Boiling CHF Model
    n_pow = 2.0d0 - 0.5d0 * (P_bar / 100.0d0)
    A_param = (2.317d0 * (D_m**0.2d0) * (hfg_kJkg) * (G_massflux**0.5d0)) / (1.0d0 + 0.0183d0 * (D_m**0.5d0) * G_massflux)
    C_param = 0.077d0 * (P_bar**0.3d0) * (D_m**0.5d0) * G_massflux

    q_CHF_kWm2 = ((A_param + 0.25d0 * D_m * G_massflux * Delta_h_sub_kJkg) / (C_param + L_heated_m)) * 0.5d0
    if (q_CHF_kWm2 < 100.0d0) q_CHF_kWm2 = 100.0d0

    DNBR_margin = q_CHF_kWm2 / max(1.0d0, q_actual_kWm2)
    Q_crit_kW = q_CHF_kWm2 * A_heat_m2

    ! Thermodynamic critical quality x_crit
    x_crit_exit = (Q_crit_kW - G_massflux * A_flow_m2 * Delta_h_sub_kJkg) / (G_massflux * A_flow_m2 * hfg_kJkg)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — FLOW BOILING CRITICAL HEAT FLUX (CHF/DNB)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Tube Diameter / Length    = ', D_tube_mm, ' mm / ', L_heated_m, ' m'
    write(*,'(A,F10.1,A,F10.1,A)') 'Mass Flux G / Pressure P  = ', G_massflux, ' kg/(m2.s) / ', P_bar, ' bar'
    write(*,'(A,F10.1,A,F10.1,A)') 'Inlet Subcool / Tsat      = ', DeltaT_sub_C, ' deg C / ', Tsat_C, ' deg C'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'CRITICAL HEAT FLUX (q_CHF)= ', q_CHF_kWm2, ' kW/m2'
    write(*,'(A,F10.2,A)')  'Applied Actual Heat Flux  = ', q_actual_kWm2, ' kW/m2'
    write(*,'(A,F10.3)')           'SAFETY MARGIN (MDNBR)     = ', DNBR_margin
    write(*,'(A,F10.2,A)')  'Critical Channel Power    = ', Q_crit_kW, ' kW'
    write(*,'(A,F10.3)')           'Critical Exit Quality x   = ', x_crit_exit
    write(*,'(A)') '============================================================'

end program flow_boiling_critical_heat_flux