program humid_air_turbine_hat_cycle
    implicit none
    integer :: iostat_val
    double precision :: r_p_overall, T_tit_C, T_amb_C, m_dot_air_dry_kgs, water_air_ratio_pct
    double precision :: eta_comp, eta_turb, T_water_in_C, eps_recup
    double precision :: m_dot_water_kgs, m_dot_humid_kgs, W_lpc_MW, W_hpc_MW, W_turb_MW
    double precision :: P_net_MW, Q_fuel_MW, eta_hat_pct, T_exhaust_C, Q_recup_MW

    ! Read inputs
    read(*,*,iostat=iostat_val) r_p_overall        ! Overall Compressor Pressure Ratio (e.g. 20.0)
    read(*,*,iostat=iostat_val) T_tit_C            ! Turbine Inlet Temp TIT [deg C] (e.g. 1200.0)
    read(*,*,iostat=iostat_val) T_amb_C            ! Ambient Air Temp [deg C] (e.g. 15.0)
    read(*,*,iostat=iostat_val) m_dot_air_dry_kgs  ! Dry Air Mass Flow [kg/s] (e.g. 100.0)
    read(*,*,iostat=iostat_val) water_air_ratio_pct! Evaporated Water Ratio mw/ma [pct] (e.g. 18.0)
    read(*,*,iostat=iostat_val) eta_comp           ! Compressor Isentropic Efficiency (e.g. 0.88)
    read(*,*,iostat=iostat_val) eta_turb           ! Turbine Isentropic Efficiency (e.g. 0.91)
    read(*,*,iostat=iostat_val) eps_recup          ! Air Recuperator Effectiveness (e.g. 0.88)

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

    if (r_p_overall <= 1.0d0 .or. m_dot_air_dry_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_water_kgs = m_dot_air_dry_kgs * (water_air_ratio_pct / 100.0d0)
    m_dot_humid_kgs = m_dot_air_dry_kgs + m_dot_water_kgs

    ! Split compression with intercooling (LPC ratio = sqrt(rp), HPC ratio = sqrt(rp))
    W_lpc_MW = (m_dot_air_dry_kgs * 1.005d0 * (T_amb_C + 273.15d0) * (((sqrt(r_p_overall))**0.286d0) - 1.0d0) / eta_comp) / 1000.0d0
    W_hpc_MW = (m_dot_air_dry_kgs * 1.005d0 * (35.0d0 + 273.15d0) * (((sqrt(r_p_overall))**0.286d0) - 1.0d0) / eta_comp) / 1000.0d0

    ! Humid Air Turbine Expansion [MW] (humid air cp ~ 1.28 kJ/kg.K)
    W_turb_MW = (m_dot_humid_kgs * 1.28d0 * (T_tit_C + 273.15d0) * (1.0d0 - (1.0d0 / (r_p_overall**0.245d0))) * eta_turb) / 1000.0d0
    P_net_MW = W_turb_MW - W_lpc_MW - W_hpc_MW

    ! Fuel Heat Duty (after recuperator preheat to ~550 C)
    Q_fuel_MW = (m_dot_humid_kgs * 1.24d0 * (T_tit_C - 550.0d0)) / 1000.0d0
    if (Q_fuel_MW < 10.0d0) Q_fuel_MW = 10.0d0
    eta_hat_pct = (P_net_MW / Q_fuel_MW) * 100.0d0

    T_exhaust_C = (T_tit_C + 273.15d0) * (1.0d0 - 0.78d0 * (1.0d0 - 1.0d0 / (r_p_overall**0.245d0))) - 273.15d0
    Q_recup_MW = (m_dot_humid_kgs * 1.20d0 * (550.0d0 - 140.0d0) * eps_recup) / 1000.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — HUMID AIR TURBINE CYCLE (HAT / EVAPORATIVE)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Pressure Ratio / TIT Temp  = ', r_p_overall, ' / ', T_tit_C, ' deg C'
    write(*,'(A,F10.1,A,F10.1,A)') 'Dry Air Flow / Humid Flow  = ', m_dot_air_dry_kgs, ' kg/s / ', m_dot_humid_kgs, ' kg/s'
    write(*,'(A,F10.2,A,F10.1,A)') 'Saturator Evaporated Water = ', m_dot_water_kgs, ' kg/s (', water_air_ratio_pct, ' %)'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'NET ELECTRICAL POWER (Pnet)= ', P_net_MW, ' MW'
    write(*,'(A,F10.2,A)')  'HAT THERMAL EFFICIENCY     = ', eta_hat_pct, ' %'
    write(*,'(A,F10.2,A)')  'Gross Turbine Expansion    = ', W_turb_MW, ' MW'
    write(*,'(A,F10.2,A)')  'Intercooled Compression Work= ', W_lpc_MW + W_hpc_MW, ' MW'
    write(*,'(A,F10.2,A)')  'Recuperator Heat Duty      = ', Q_recup_MW, ' MW'
    write(*,'(A,F10.1,A)')  'Turbine Exhaust Gas Temp   = ', T_exhaust_C, ' deg C'
    write(*,'(A)') '============================================================'

end program humid_air_turbine_hat_cycle