Supercritical CO2 Recompression Brayton

program supercritical_co2_brayton_recompression
    implicit none
    integer :: iostat_val
    double precision :: P_high_MPa, P_low_MPa, T_tit_C, T_mc_in_C, m_dot_total_kgs, recomp_frac
    double precision :: eta_turb, eta_mc, eta_rc, eps_ltr, eps_htr
    double precision :: w_turb, w_mc, w_rc, w_net, q_heater, eta_th, P_net_MW, Q_heater_MW, Q_precooler_MW
    double precision :: T_mc_out_C, T_rc_out_C, T_turb_out_C

    ! Read inputs
    read(*,*,iostat=iostat_val) P_high_MPa       ! High Pressure [MPa] (e.g. 25.0)
    read(*,*,iostat=iostat_val) P_low_MPa        ! Low Pressure (MC Inlet) [MPa] (e.g. 7.8)
    read(*,*,iostat=iostat_val) T_tit_C          ! Turbine Inlet Temp [deg C] (e.g. 650.0)
    read(*,*,iostat=iostat_val) T_mc_in_C        ! Main Compressor Inlet Temp [deg C] (e.g. 32.0)
    read(*,*,iostat=iostat_val) m_dot_total_kgs  ! Total CO2 Mass Flow [kg/s] (e.g. 100.0)
    read(*,*,iostat=iostat_val) recomp_frac      ! Recompression Fraction gamma (e.g. 0.32)
    read(*,*,iostat=iostat_val) eta_turb         ! Turbine Efficiency (e.g. 0.90)
    read(*,*,iostat=iostat_val) eta_mc           ! Main Compressor Efficiency (e.g. 0.86)
    read(*,*,iostat=iostat_val) eta_rc           ! Recompressor Efficiency (e.g. 0.86)

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

    if (P_high_MPa <= P_low_MPa .or. T_tit_C <= T_mc_in_C .or. m_dot_total_kgs <= 0.0d0) then
        write(*,*) 'ERROR: P_high > P_low and T_tit > T_mc_in required.'
        stop
    end if

    ! Enthalpy differences & Specific Work (kJ/kg) for sCO2
    ! Turbine expansion from P_high to P_low
    w_turb = eta_turb * (1.22d0 * (T_tit_C + 273.15d0) * (1.0d0 - (P_low_MPa / P_high_MPa)**0.24d0))

    ! Main Compressor (near-critical liquid-like dense CO2: rho ~ 600 kg/m3)
    w_mc = (((P_high_MPa - P_low_MPa) * 1000.0d0 / 580.0d0) + 12.0d0) / eta_mc

    ! Recompressor (vapor-like CO2 from LTR outlet)
    w_rc = (1.15d0 * 350.0d0 * ((P_high_MPa / P_low_MPa)**0.22d0 - 1.0d0)) / eta_rc

    ! Heat addition in main heater (after HTR preheating)
    q_heater = 1.25d0 * (T_tit_C - 480.0d0) + 180.0d0
    if (q_heater < 200.0d0) q_heater = 200.0d0

    ! Net Cycle Work
    w_net = w_turb - ((1.0d0 - recomp_frac) * w_mc) - (recomp_frac * w_rc)
    eta_th = (w_net / q_heater) * 100.0d0

    P_net_MW = (m_dot_total_kgs * w_net) / 1000.0d0
    Q_heater_MW = (m_dot_total_kgs * q_heater) / 1000.0d0
    Q_precooler_MW = Q_heater_MW - P_net_MW

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERCRITICAL CO2 RECOMPRESSION BRAYTON'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'P_high / P_low (MC Inlet) = ', P_high_MPa, ' MPa / ', P_low_MPa, ' MPa'
    write(*,'(A,F10.1,A,F10.1,A)') 'TIT Temp / MC Inlet Temp  = ', T_tit_C, ' deg C / ', T_mc_in_C, ' deg C'
    write(*,'(A,F10.2,A,F10.2)')   'Mass Flow / Recomp Frac   = ', m_dot_total_kgs, ' kg/s / ', recomp_frac
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'NET ELECTRICAL POWER      = ', P_net_MW, ' MW'
    write(*,'(A,F10.2,A)')  'THERMAL EFFICIENCY (eta)  = ', eta_th, ' %'
    write(*,'(A,F10.2,A)')  'Turbine Power Output      = ', (m_dot_total_kgs*w_turb)/1000.0d0, ' MW'
    write(*,'(A,F10.2,A)')  'Main + Recomp Comp Power  = ', (m_dot_total_kgs*((1.0d0-recomp_frac)*w_mc + recomp_frac*w_rc))/1000.0d0, ' MW'
    write(*,'(A,F10.2,A)')  'Reactor / Heater Duty     = ', Q_heater_MW, ' MW'
    write(*,'(A,F10.2,A)')  'Precooler Rejection Duty  = ', Q_precooler_MW, ' MW'
    write(*,'(A)') '============================================================'

end program supercritical_co2_brayton_recompression