Transcritical CO2 Refrigeration (R744)

program transcritical_co2_refrigeration_r744
    implicit none
    integer :: iostat_val
    double precision :: Te_evap_C, T_gcout_C, P_high_bar, m_dot_kgs, eta_comp, DT_superheat_K
    double precision :: Pe_bar, P_opt_bar, h1, h2, h3, h4, q_evap, w_comp, cop_cool, cop_heat
    double precision :: Q_cool_kW, W_elec_kW, Q_gas_cooler_kW, T_comp_disch_C

    ! Read inputs
    read(*,*,iostat=iostat_val) Te_evap_C       ! Evaporating Temp [deg C] (e.g. -8.0 for Medium Temp MT or -30.0 for LT)
    read(*,*,iostat=iostat_val) T_gcout_C      ! Gas Cooler Outlet Temp [deg C] (e.g. 38.0)
    read(*,*,iostat=iostat_val) P_high_bar     ! High-Side Supercritical Pressure [bar] (0 for auto-optimum) (e.g. 95.0)
    read(*,*,iostat=iostat_val) m_dot_kgs      ! R744 Mass Flow [kg/s] (e.g. 1.20)
    read(*,*,iostat=iostat_val) eta_comp       ! Compressor Isentropic Efficiency (e.g. 0.72)
    read(*,*,iostat=iostat_val) DT_superheat_K ! Suction Superheat [K] (e.g. 5.0)

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

    ! Approximate R744 saturation pressure at Te: Pe [bar]
    Pe_bar = 30.4d0 + 0.95d0 * Te_evap_C + 0.007d0 * (Te_evap_C**2)
    if (Pe_bar < 10.0d0) Pe_bar = 10.0d0

    ! Optimum High Pressure (Liao / Kauf correlation)
    P_opt_bar = 2.6d0 * T_gcout_C + 7.5d0
    if (P_opt_bar < 75.0d0) P_opt_bar = 75.0d0

    if (P_high_bar <= 74.0d0) then
        P_high_bar = P_opt_bar
    end if

    ! Enthalpies (kJ/kg) for R744 near transcritical dome
    h1 = 430.0d0 + 0.85d0 * Te_evap_C + 1.8d0 * DT_superheat_K
    ! Gas cooler exit enthalpy (strong function of P_high and T_gcout)
    h3 = 190.0d0 + 3.2d0 * (T_gcout_C - 25.0d0) - 0.75d0 * (P_high_bar - 80.0d0)
    if (h3 < 150.0d0) h3 = 150.0d0
    h4 = h3 ! isenthalpic expansion valve

    ! Compression work
    w_comp = ((1.25d0 * (P_high_bar / Pe_bar)**0.286d0 - 1.0d0) * 160.0d0) / eta_comp
    h2 = h1 + w_comp
    T_comp_disch_C = Te_evap_C + DT_superheat_K + (w_comp / 1.35d0)

    ! Performance metrics
    q_evap = max(10.0d0, h1 - h4)
    cop_cool = q_evap / w_comp
    cop_heat = (h2 - h3) / w_comp

    Q_cool_kW = m_dot_kgs * q_evap
    W_elec_kW = m_dot_kgs * w_comp
    Q_gas_cooler_kW = m_dot_kgs * (h2 - h3)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TRANSCRITICAL CO2 (R744) REFRIGERATION'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Evap Temp Te / Sat Press Pe= ', Te_evap_C, ' deg C / ', Pe_bar, ' bar'
    write(*,'(A,F10.2,A,F10.2,A)') 'Gas Cooler Tout / High P   = ', T_gcout_C, ' deg C / ', P_high_bar, ' bar'
    write(*,'(A,F10.2,A)')         'Optimal High-Side Pressure = ', P_opt_bar, ' bar'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'COOLING CAPACITY (Q_cool) = ', Q_cool_kW, ' kW'
    write(*,'(A,F10.2,A)')  'COMPRESSOR POWER (W_elec) = ', W_elec_kW, ' kW'
    write(*,'(A,F10.3)')    'COOLING COP (COP_cool)    = ', cop_cool
    write(*,'(A,F10.3)')    'HEATING COP (COP_heat)    = ', cop_heat
    write(*,'(A,F10.2,A)')  'Gas Cooler Rejection Heat = ', Q_gas_cooler_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Compressor Discharge Temp = ', T_comp_disch_C, ' deg C'
    write(*,'(A)') '============================================================'

end program transcritical_co2_refrigeration_r744