program microchannel_flow_boiling
    implicit none
    integer :: iostat_val, N_channels
    double precision :: W_ch_mm, H_ch_mm, L_ch_mm, G_massflux, q_flux_Wcm2, x_quality, Psat_bar
    double precision :: rho_L, rho_V, mu_L_cP, k_L, cp_L, hfg_kJkg, sigma_Nm, F_fluid
    double precision :: W_m, H_m, L_m, Dh_m, Ac_m2, mu_L_Pas, Co_num, Re_LO, Pr_L, h_LO
    double precision :: Bo_num, q_flux_Wm2, h_NBD, h_CBD, h_tp_Wm2K, Tsat_C, Twall_C, DeltaP_tp_kPa
    double precision, parameter :: G_GRAV = 9.80665d0

    ! Read inputs
    read(*,*,iostat=iostat_val) W_ch_mm        ! Channel Width [mm] (e.g. 0.30)
    read(*,*,iostat=iostat_val) H_ch_mm        ! Channel Height [mm] (e.g. 0.80)
    read(*,*,iostat=iostat_val) L_ch_mm        ! Channel Length [mm] (e.g. 20.0)
    read(*,*,iostat=iostat_val) N_channels     ! Number of Parallel Microchannels (e.g. 50)
    read(*,*,iostat=iostat_val) G_massflux     ! Mass Flux G [kg/(m2.s)] (e.g. 600.0)
    read(*,*,iostat=iostat_val) q_flux_Wcm2    ! Heat Flux [W/cm2] (e.g. 80.0)
    read(*,*,iostat=iostat_val) x_quality      ! Mean Vapor Quality (0.05 to 0.9) (e.g. 0.35)
    read(*,*,iostat=iostat_val) Psat_bar       ! Saturation Pressure [bar] (e.g. 7.5 for R134a, 1.0 for Water)
    read(*,*,iostat=iostat_val) rho_L          ! Liquid Density [kg/m3] (e.g. 1200.0)
    read(*,*,iostat=iostat_val) rho_V          ! Vapor Density [kg/m3] (e.g. 35.0)
    read(*,*,iostat=iostat_val) mu_L_cP        ! Liquid Viscosity [cP] (e.g. 0.19)
    read(*,*,iostat=iostat_val) k_L            ! Liquid Conductivity [W/(m.K)] (e.g. 0.082)
    read(*,*,iostat=iostat_val) cp_L           ! Liquid Specific Heat [J/(kg.K)] (e.g. 1430.0)
    read(*,*,iostat=iostat_val) hfg_kJkg       ! Latent Heat [kJ/kg] (e.g. 175.0)
    read(*,*,iostat=iostat_val) sigma_Nm       ! Surface Tension [N/m] (e.g. 0.008)
    read(*,*,iostat=iostat_val) F_fluid        ! Kandlikar Fluid Factor (e.g. 1.63 for R134a, 1.0 for Water)

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

    if (W_ch_mm <= 0.0d0 .or. H_ch_mm <= 0.0d0 .or. G_massflux <= 0.0d0 .or. q_flux_Wcm2 <= 0.0d0) then
        write(*,*) 'ERROR: Dimensions, mass flux, and heat flux must be positive.'
        stop
    end if

    W_m = W_ch_mm * 1.0d-3
    H_m = H_ch_mm * 1.0d-3
    L_m = L_ch_mm * 1.0d-3
    Dh_m = (2.0d0 * W_m * H_m) / (W_m + H_m)
    Ac_m2 = W_m * H_m
    mu_L_Pas = mu_L_cP * 1.0d-3
    q_flux_Wm2 = q_flux_Wcm2 * 10000.0d0 ! W/m2

    ! Confinement Number
    Co_num = (1.0d0 / Dh_m) * sqrt(sigma_Nm / (G_GRAV * (rho_L - rho_V)))

    ! Liquid-only Reynolds and Nusselt
    Re_LO = (G_massflux * Dh_m) / mu_L_Pas
    Pr_L = (cp_L * mu_L_Pas) / k_L
    h_LO = (0.023d0 * (Re_LO**0.8d0) * (Pr_L**0.4d0) * k_L) / Dh_m

    ! Boiling Number
    Bo_num = q_flux_Wm2 / (G_massflux * hfg_kJkg * 1000.0d0)

    ! Kandlikar (2004) Microchannel Flow Boiling HTC
    h_NBD = 0.6683d0 * (max(0.1d0, Co_num)**(-0.2d0)) * ((1.0d0 - x_quality)**0.8d0) * h_LO + &
            1058.0d0 * (Bo_num**0.7d0) * ((1.0d0 - x_quality)**0.8d0) * F_fluid * h_LO

    h_CBD = 1.136d0 * (max(0.1d0, Co_num)**(-0.9d0)) * ((1.0d0 - x_quality)**0.8d0) * h_LO + &
            667.0d0 * (Bo_num**0.7d0) * ((1.0d0 - x_quality)**0.8d0) * F_fluid * h_LO

    h_tp_Wm2K = max(h_NBD, h_CBD)
    if (h_tp_Wm2K < 500.0d0) h_tp_Wm2K = 500.0d0

    Tsat_C = 28.0d0 * (Psat_bar**0.35d0) ! approx saturation temp
    Twall_C = Tsat_C + (q_flux_Wm2 / h_tp_Wm2K)

    ! Two-Phase Pressure Drop DeltaP [kPa]
    DeltaP_tp_kPa = ((0.079d0 / (Re_LO**0.25d0)) * (L_m / Dh_m) * (G_massflux**2) / (2.0d0 * rho_L) * &
                    (1.0d0 + x_quality * (rho_L / rho_V - 1.0d0))) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — MICROCHANNEL FLOW BOILING ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.3,A,I5,A)')    'Hydraulic Diam Dh / N_ch   = ', Dh_m*1000.0d0, ' mm / ', N_channels, ' channels'
    write(*,'(A,F10.3,A,E12.4)')   'Confinement Co / Boiling Bo= ', Co_num, ' / ', Bo_num
    write(*,'(A,F10.1,A,F10.1,A)') 'Mass Flux G / Heat Flux q" = ', G_massflux, ' kg/(m2.s) / ', q_flux_Wcm2, ' W/cm2'
    write(*,'(A,F10.2,A,F10.2,A)') 'Vapor Quality x / Psat     = ', x_quality, ' / ', Psat_bar, ' bar'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.1,A)')  'TWO-PHASE BOILING HTC h_tp= ', h_tp_Wm2K, ' W/(m2.K)'
    write(*,'(A,F10.1,A)')  'Liquid-Only HTC h_LO      = ', h_LO, ' W/(m2.K)'
    write(*,'(A,F10.2,A,F10.2,A)') 'Channel Wall / Saturation T= ', Twall_C, ' deg C / ', Tsat_C, ' deg C'
    write(*,'(A,F10.2,A)')  'Two-Phase Pressure Drop   = ', DeltaP_tp_kPa, ' kPa'
    write(*,'(A)') '============================================================'

end program microchannel_flow_boiling