program ejector_refrigeration_solar_thermal
    implicit none
    integer :: iostat_val, fluid_type
    double precision :: Tg_gen_C, Te_evap_C, Tc_cond_C, Q_cool_kW, I_solar_Wm2, eta_collector
    double precision :: Pg_bar, Pe_bar, Pc_bar, entrain_ratio_omega, cop_thermal, Q_gen_kW, A_solar_m2
    double precision :: m_dot_pri_kgs, m_dot_sec_kgs, W_pump_W

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type       ! 1=Water(R718), 2=R134a, 3=R245fa, 4=Isobutane(R600a)
    read(*,*,iostat=iostat_val) Tg_gen_C         ! Solar Generator/Boiler Temp [deg C] (e.g. 90.0)
    read(*,*,iostat=iostat_val) Te_evap_C        ! Evaporator Temp [deg C] (e.g. 8.0)
    read(*,*,iostat=iostat_val) Tc_cond_C        ! Condenser Temp [deg C] (e.g. 35.0)
    read(*,*,iostat=iostat_val) Q_cool_kW        ! Target Cooling Capacity [kW] (e.g. 10.0)
    read(*,*,iostat=iostat_val) I_solar_Wm2      ! Solar Global Irradiance [W/m2] (e.g. 850.0)
    read(*,*,iostat=iostat_val) eta_collector    ! Solar Collector Thermal Efficiency (e.g. 0.60)

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

    if (Tg_gen_C <= Tc_cond_C .or. Tc_cond_C <= Te_evap_C .or. Q_cool_kW <= 0.0d0) then
        write(*,*) 'ERROR: Requires Tg > Tc > Te and Q_cool > 0.'
        stop
    end if

    ! Approximate Saturation Pressures
    if (fluid_type == 1) then ! Water R718
        Pg_bar = 0.70d0 * (Tg_gen_C / 90.0d0)**3.5d0
        Pe_bar = 0.0107d0 ! at 8 deg C
        Pc_bar = 0.0562d0 ! at 35 deg C
    else if (fluid_type == 2) then ! R134a
        Pg_bar = 32.0d0 * (Tg_gen_C / 90.0d0)**2.5d0
        Pe_bar = 3.88d0
        Pc_bar = 8.87d0
    else if (fluid_type == 3) then ! R245fa
        Pg_bar = 10.2d0 * (Tg_gen_C / 90.0d0)**2.8d0
        Pe_bar = 0.76d0
        Pc_bar = 2.12d0
    else ! Isobutane
        Pg_bar = 15.0d0 * (Tg_gen_C / 90.0d0)**2.6d0
        Pe_bar = 2.05d0
        Pc_bar = 4.70d0
    end if

    ! Huang (1999) 1D Supersonic Ejector Entrainment Ratio Model
    entrain_ratio_omega = 0.55d0 * sqrt((Tg_gen_C - Tc_cond_C) / (Tc_cond_C - Te_evap_C + 0.1d0)) * &
                          ((Pe_bar / Pc_bar)**0.45d0)
    if (entrain_ratio_omega < 0.10d0) entrain_ratio_omega = 0.10d0
    if (entrain_ratio_omega > 1.20d0) entrain_ratio_omega = 1.20d0

    ! Thermal COP of Ejector Cycle
    cop_thermal = entrain_ratio_omega * 0.85d0 * ((Te_evap_C + 273.15d0) / (Tg_gen_C - Te_evap_C + 10.0d0))
    if (cop_thermal > 0.85d0) cop_thermal = 0.85d0
    if (cop_thermal < 0.15d0) cop_thermal = 0.15d0

    ! Solar Thermal Heat Requirements
    Q_gen_kW = Q_cool_kW / cop_thermal
    A_solar_m2 = (Q_gen_kW * 1000.0d0) / (max(100.0d0, I_solar_Wm2) * max(0.2d0, eta_collector))

    m_dot_sec_kgs = Q_cool_kW / 220.0d0 ! approx latent heat
    m_dot_pri_kgs = m_dot_sec_kgs / entrain_ratio_omega
    W_pump_W = m_dot_pri_kgs * (Pg_bar - Pc_bar) * 100.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SOLAR EJECTOR REFRIGERATION CYCLE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Generator Tg / Evap Te     = ', Tg_gen_C, ' deg C / ', Te_evap_C, ' deg C'
    write(*,'(A,F10.1,A,F10.2,A)') 'Condenser Tc / Cooling Q_c = ', Tc_cond_C, ' deg C / ', Q_cool_kW, ' kW'
    write(*,'(A,F10.3,A,F10.3)')   'Entrainment Ratio omega / COP= ', entrain_ratio_omega, ' / ', cop_thermal
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'SOLAR COLLECTOR AREA      = ', A_solar_m2, ' m2'
    write(*,'(A,F10.2,A)')  'Solar Generator Heat Duty = ', Q_gen_kW, ' kW'
    write(*,'(A,F10.3,A,F10.3,A)') 'Motive Flow / Suction Flow= ', m_dot_pri_kgs, ' kg/s / ', m_dot_sec_kgs, ' kg/s'
    write(*,'(A,F10.1,A)')  'Feed Pump Power (Electric)= ', W_pump_W, ' Watts'
    write(*,'(A)') '============================================================'

end program ejector_refrigeration_solar_thermal