Turbocharger Engine Match & Surge Margin Sizer

program turbocharger_match_engine_surge_margin
    implicit none
    integer :: iostat_val, engine_type
    double precision :: Vd_liters, N_engine_rpm, P_boost_bar_g, T_man_C, eta_vol_eng
    double precision :: Texh_C, eta_comp, eta_turb, eta_mech
    double precision :: P1_bar, P_man_bar, T1_K, Tman_K, Texh_K, m_dot_air_kgs, m_dot_exh_kgs
    double precision :: PR_c, P_comp_kW, P_turb_req_kW, PR_t, P_exh_bar_g
    double precision :: m_dot_corr_kgs, surge_mass_flow, surge_margin_pct, N_turbo_rpm
    character(len=32) :: match_status
    double precision, parameter :: R_AIR = 287.05d0
    double precision, parameter :: CP_AIR = 1005.0d0
    double precision, parameter :: GAMMA_AIR = 1.40d0
    double precision, parameter :: CP_EXH = 1150.0d0
    double precision, parameter :: GAMMA_EXH = 1.33d0

    ! Read inputs
    read(*,*,iostat=iostat_val) engine_type      ! 1=Gasoline Turbo, 2=Diesel HD Truck, 3=Marine Propulsion, 4=Motorsport
    read(*,*,iostat=iostat_val) Vd_liters        ! Engine Displacement [liters] (e.g. 2.0)
    read(*,*,iostat=iostat_val) N_engine_rpm     ! Engine Rated RPM [rpm] (e.g. 6000.0)
    read(*,*,iostat=iostat_val) P_boost_bar_g    ! Boost Pressure [bar gauge] (e.g. 1.40)
    read(*,*,iostat=iostat_val) T_man_C          ! Intercooled Manifold Temp [deg C] (e.g. 45.0)
    read(*,*,iostat=iostat_val) eta_vol_eng      ! Engine Volumetric Efficiency [0.80 to 0.98] (e.g. 0.92)
    read(*,*,iostat=iostat_val) Texh_C           ! Turbine Inlet Exhaust Temp [deg C] (e.g. 900.0)
    read(*,*,iostat=iostat_val) eta_comp         ! Compressor Efficiency [0.65 to 0.82] (e.g. 0.74)
    read(*,*,iostat=iostat_val) eta_turb         ! Turbine Efficiency [0.65 to 0.82] (e.g. 0.72)

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid input data for turbocharger engine match.'
        stop
    end if

    if (Vd_liters <= 0.0d0 .or. N_engine_rpm <= 0.0d0 .or. P_boost_bar_g <= 0.0d0) then
        write(*,*) 'ERROR: Engine displacement, RPM, and boost must be positive.'
        stop
    end if

    P1_bar = 1.013d0
    P_man_bar = P1_bar + P_boost_bar_g
    T1_K = 293.15d0
    Tman_K = T_man_C + 273.15d0
    Texh_K = Texh_C + 273.15d0

    ! Engine airflow (4-stroke) [kg/s]
    m_dot_air_kgs = ((Vd_liters * 1.0d-3) * (N_engine_rpm / 120.0d0) * eta_vol_eng * (P_man_bar * 1.0d5)) / (R_AIR * Tman_K)
    m_dot_exh_kgs = m_dot_air_kgs * 1.065d0 ! fuel addition (A/F ~ 15)

    PR_c = P_man_bar / (P1_bar - 0.03d0) ! 30 mbar airbox loss

    ! Compressor Power Demand [kW]
    P_comp_kW = (m_dot_air_kgs * CP_AIR * T1_K / 1000.0d0 / max(0.50d0, eta_comp)) * &
                (PR_c**((GAMMA_AIR - 1.0d0) / GAMMA_AIR) - 1.0d0)

    ! Mechanical efficiency of bearings (journal/ball)
    eta_mech = 0.95d0
    P_turb_req_kW = P_comp_kW / eta_mech

    ! Turbine Expansion Ratio & Backpressure
    PR_t = (1.0d0 / max(0.20d0, 1.0d0 - (P_turb_req_kW * 1000.0d0) / &
           (m_dot_exh_kgs * CP_EXH * Texh_K * max(0.50d0, eta_turb))))**(GAMMA_EXH / (GAMMA_EXH - 1.0d0))
    PR_t = min(4.5d0, max(1.15d0, PR_t))

    P_exh_bar_g = (P1_bar * PR_t) - P1_bar

    ! Corrected Mass Flow Rate
    m_dot_corr_kgs = m_dot_air_kgs * sqrt(T1_K / 288.15d0) / (P1_bar / 1.01325d0)

    ! Surge Line Modeling (empirical surge flow at PR_c)
    surge_mass_flow = 0.35d0 * m_dot_corr_kgs * (PR_c / 2.0d0)**0.85d0
    surge_margin_pct = ((m_dot_corr_kgs - surge_mass_flow) / max(0.01d0, m_dot_corr_kgs)) * 100.0d0
    surge_margin_pct = max(5.0d0, min(45.0d0, surge_margin_pct))

    ! Turbocharger Shaft Speed (Estimated from tip diameter ~45-75mm)
    N_turbo_rpm = 165000.0d0 * sqrt(PR_c / 2.2d0) * (2.0d0 / max(0.5d0, Vd_liters))**0.25d0
    N_turbo_rpm = min(260000.0d0, max(60000.0d0, N_turbo_rpm))

    if (surge_margin_pct < 12.0d0) then
        match_status = 'SURGE DANGER ZONE (SM < 12%)'
    else if (surge_margin_pct <= 28.0d0) then
        match_status = 'OPTIMAL MATCHED CRUISE'
    else
        match_status = 'CHOKE / HIGH FLOW REGIME'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — TURBOCHARGER ENGINE MATCH & SURGE SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Compressor PR / Boost        = ', PR_c, ' / ', P_boost_bar_g, ' bar g'
    write(*,'(A,F10.2,A,F10.1,A)') 'Turbine PR / Exhaust Backpr. = ', PR_t, ' / ', P_exh_bar_g, ' bar g'
    write(*,'(A,F10.1,A,A)')       'Surge Margin                 = ', surge_margin_pct, ' % | ', trim(match_status)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A,F10.0,A)') 'COMPRESSOR POWER DEMAND      = ', P_comp_kW, ' kW (Rotor: ', N_turbo_rpm, ' rpm)'
    write(*,'(A,F10.3,A,F10.1,A)') 'Engine Air Mass Flow         = ', m_dot_air_kgs, ' kg/s (', m_dot_air_kgs*3600.0d0, ' kg/h)'
    write(*,'(A,F10.1,A,F10.1,A)') 'Exhaust Temp / Manifold Temp = ', Texh_C, ' C / ', T_man_C, ' C'
    write(*,'(A,F10.2,A,F10.2,A)') 'Compressor / Turbine Effic.  = ', eta_comp*100.0d0, ' % / ', eta_turb*100.0d0, ' %'
    write(*,'(A)') '============================================================'

end program turbocharger_match_engine_surge_margin