Steam-Injected Gas Turbine (STIG / Cheng)

program steam_injected_gas_turbine_stig
    implicit none
    integer :: iostat_val
    double precision :: r_p_ratio, T_tit_C, T_amb_C, m_dot_air_kgs, steam_ratio_pct
    double precision :: eta_comp, eta_turb, eta_comb, T_steam_inj_C, P_inj_bar
    double precision :: m_dot_steam_kgs, W_comp_MW, W_turb_simple_MW, W_turb_stig_MW
    double precision :: P_net_simple_MW, P_net_stig_MW, eta_simple_pct, eta_stig_pct
    double precision :: Q_fuel_MW, power_boost_pct, T_exhaust_C, Q_hrsg_MW

    ! Read inputs
    read(*,*,iostat=iostat_val) r_p_ratio        ! Compressor Pressure Ratio (e.g. 16.0)
    read(*,*,iostat=iostat_val) T_tit_C          ! Turbine Inlet Temp TIT [deg C] (e.g. 1150.0)
    read(*,*,iostat=iostat_val) T_amb_C          ! Ambient Air Temp [deg C] (e.g. 20.0)
    read(*,*,iostat=iostat_val) m_dot_air_kgs    ! Inlet Air Mass Flow [kg/s] (e.g. 80.0)
    read(*,*,iostat=iostat_val) steam_ratio_pct  ! Steam Injection Ratio ms/ma [pct] (e.g. 10.0)
    read(*,*,iostat=iostat_val) eta_comp         ! Compressor Isentropic Efficiency (e.g. 0.86)
    read(*,*,iostat=iostat_val) eta_turb         ! Turbine Isentropic Efficiency (e.g. 0.89)
    read(*,*,iostat=iostat_val) T_steam_inj_C    ! Superheated Steam Temp [deg C] (e.g. 450.0)

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

    if (r_p_ratio <= 1.0d0 .or. m_dot_air_kgs <= 0.0d0 .or. T_tit_C <= T_amb_C) then
        write(*,*) 'ERROR: Pressure ratio > 1 and TIT > Tamb required.'
        stop
    end if

    m_dot_steam_kgs = m_dot_air_kgs * (steam_ratio_pct / 100.0d0)
    P_inj_bar = r_p_ratio * 1.15d0

    ! Compressor Work (Air only) [MW]
    W_comp_MW = (m_dot_air_kgs * 1.005d0 * (T_amb_C + 273.15d0) * ((r_p_ratio**0.286d0) - 1.0d0) / eta_comp) / 1000.0d0

    ! Simple Cycle Turbine Work [MW]
    W_turb_simple_MW = (m_dot_air_kgs * 1.15d0 * (T_tit_C + 273.15d0) * (1.0d0 - (1.0d0 / (r_p_ratio**0.24d0))) * eta_turb) / 1000.0d0
    P_net_simple_MW = W_turb_simple_MW - W_comp_MW

    ! Fuel Heat Duty in Combustor [MW]
    Q_fuel_MW = (m_dot_air_kgs * 1.10d0 * (T_tit_C - (T_amb_C + 350.0d0))) / 1000.0d0
    if (Q_fuel_MW < 10.0d0) Q_fuel_MW = 10.0d0
    eta_simple_pct = (P_net_simple_MW / Q_fuel_MW) * 100.0d0

    ! STIG Dual-Fluid Turbine Work (higher mass flow + higher mixture cp ~ 1.35 kJ/kg.K) [MW]
    W_turb_stig_MW = ((m_dot_air_kgs * 1.15d0 + m_dot_steam_kgs * 2.15d0) * (T_tit_C + 273.15d0) * &
                     (1.0d0 - (1.0d0 / (r_p_ratio**0.24d0))) * eta_turb) / 1000.0d0
    P_net_stig_MW = W_turb_stig_MW - W_comp_MW

    ! Additional fuel to heat injected steam to TIT
    Q_fuel_MW = Q_fuel_MW + (m_dot_steam_kgs * 2.10d0 * (T_tit_C - T_steam_inj_C)) / 1000.0d0
    eta_stig_pct = (P_net_stig_MW / Q_fuel_MW) * 100.0d0

    power_boost_pct = ((P_net_stig_MW - P_net_simple_MW) / P_net_simple_MW) * 100.0d0
    T_exhaust_C = (T_tit_C + 273.15d0) * (1.0d0 - 0.75d0 * (1.0d0 - 1.0d0 / (r_p_ratio**0.24d0))) - 273.15d0
    Q_hrsg_MW = (m_dot_steam_kgs * (3200.0d0 - 200.0d0)) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — STEAM-INJECTED GAS TURBINE (STIG / CHENG)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Pressure Ratio / TIT Temp  = ', r_p_ratio, ' / ', T_tit_C, ' deg C'
    write(*,'(A,F10.1,A,F10.1,A)') 'Air Flow / Steam Injection = ', m_dot_air_kgs, ' kg/s / ', steam_ratio_pct, ' %'
    write(*,'(A,F10.2,A)')         'Steam Injection Mass Flow  = ', m_dot_steam_kgs, ' kg/s'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.2,A)') 'STIG POWER / SIMPLE POWER  = ', P_net_stig_MW, ' MW / ', P_net_simple_MW, ' MW'
    write(*,'(A,F10.2,A)')  'POWER BOOST GAIN           = +', power_boost_pct, ' %'
    write(*,'(A,F10.2,A,F10.2,A)') 'STIG ETA / SIMPLE ETA      = ', eta_stig_pct, ' % / ', eta_simple_pct, ' %'
    write(*,'(A,F10.2,A)')  'HRSG Steam Boiler Heat Duty= ', Q_hrsg_MW, ' MW'
    write(*,'(A,F10.1,A)')  'Turbine Exhaust Gas Temp   = ', T_exhaust_C, ' deg C'
    write(*,'(A)') '============================================================'

end program steam_injected_gas_turbine_stig