Steam Turbine Wilson Loss & Baumann Sizer

program multistage_steam_turbine_wilson_loss
    implicit none
    integer :: iostat_val, num_stages
    double precision :: m_dot_steam_kgs, P1_bar, T1_C, P2_bar, eta_dry_target, baumann_factor
    double precision :: h1_kJkg, s1_kJkgK, h2s_kJkg, h2_kJkg, x2_quality, y2_moisture_pct
    double precision :: h_f_kJkg, h_g_kJkg, s_f_kJkgK, s_g_kJkgK, T_sat_C
    double precision :: delta_h_kJkg, eta_wet_actual, P_turb_MW, reheat_factor
    character(len=32) :: erosion_risk
    double precision, parameter :: PI = 3.141592653589793d0

    ! Read inputs
    read(*,*,iostat=iostat_val) m_dot_steam_kgs  ! Steam Flow Rate [kg/s] (e.g. 150.0)
    read(*,*,iostat=iostat_val) P1_bar           ! Turbine Inlet Pressure P1 [bar a] (e.g. 40.0)
    read(*,*,iostat=iostat_val) T1_C             ! Turbine Inlet Temperature T1 [deg C] (e.g. 450.0)
    read(*,*,iostat=iostat_val) P2_bar           ! Condenser Backpressure P2 [bar a] (e.g. 0.05)
    read(*,*,iostat=iostat_val) num_stages       ! Number of LP Stages (e.g. 6)
    read(*,*,iostat=iostat_val) eta_dry_target   ! Base Dry Stage Efficiency [0.80 to 0.94] (e.g. 0.90)
    read(*,*,iostat=iostat_val) baumann_factor   ! Baumann Loss Factor α [0.80 to 1.10] (e.g. 0.95)

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

    if (m_dot_steam_kgs <= 0.0d0 .or. P1_bar <= 0.0d0 .or. P2_bar <= 0.0d0 .or. P2_bar >= P1_bar) then
        write(*,*) 'ERROR: Flow rate and pressures must be positive with P2 < P1.'
        stop
    end if

    ! Inlet Steam State (Approximate Superheated Steam IAPWS-IF97 fit)
    h1_kJkg = 2500.0d0 + 2.10d0 * T1_C - 0.85d0 * P1_bar
    s1_kJkgK = 6.45d0 + 0.0035d0 * T1_C - 0.35d0 * log(max(0.1d0, P1_bar))

    ! Condenser Saturated State at P2 (Approximate)
    T_sat_C = 99.6d0 * (P2_bar**0.25d0)
    h_f_kJkg = 4.18d0 * T_sat_C
    h_g_kJkg = 2501.0d0 + 1.85d0 * T_sat_C
    s_f_kJkgK = (4.18d0 * log((T_sat_C + 273.15d0) / 273.15d0))
    s_g_kJkgK = s_f_kJkgK + (h_g_kJkg - h_f_kJkg) / (T_sat_C + 273.15d0)

    ! Isentropic Expansion Quality at P2
    x2_quality = (s1_kJkgK - s_f_kJkgK) / max(0.1d0, (s_g_kJkgK - s_f_kJkgK))
    x2_quality = min(1.0d0, max(0.70d0, x2_quality))
    h2s_kJkg = h_f_kJkg + x2_quality * (h_g_kJkg - h_f_kJkg)

    ! Baumann Rule: Wet Efficiency Degradation
    y2_moisture_pct = (1.0d0 - x2_quality) * 100.0d0
    eta_wet_actual = eta_dry_target * (1.0d0 - baumann_factor * (y2_moisture_pct / 200.0d0))
    eta_wet_actual = max(0.65d0, min(eta_dry_target, eta_wet_actual))

    ! Multi-stage Reheat Factor
    reheat_factor = 1.0d0 + 0.035d0 * dble(num_stages) * (1.0d0 - eta_dry_target)
    delta_h_kJkg = (h1_kJkg - h2s_kJkg) * eta_wet_actual * reheat_factor
    h2_kJkg = h1_kJkg - delta_h_kJkg

    P_turb_MW = (m_dot_steam_kgs * delta_h_kJkg) / 1000.0d0

    ! Blade Erosion Risk Classification
    if (y2_moisture_pct < 8.0d0) then
        erosion_risk = 'VERY LOW (DRY/EARLY WET)'
    else if (y2_moisture_pct <= 12.0d0) then
        erosion_risk = 'MODERATE (STANDARD DRAINAGE)'
    else
        erosion_risk = 'HIGH (STELLITE SHIELDING REQ)'
    end if

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — STEAM TURBINE WILSON & BAUMANN SIZER'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'Inlet / Exit Pressure        = ', P1_bar, ' bar / ', P2_bar, ' bar'
    write(*,'(A,F10.1,A,F10.2,A)') 'Inlet Temp / Exit Quality x2 = ', T1_C, ' C / ', x2_quality, ''
    write(*,'(A,F10.1,A,A)')       'Exhaust Moisture (y2)        = ', y2_moisture_pct, ' % | ', trim(erosion_risk)
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'TURBINE POWER GENERATED      = ', P_turb_MW, ' MW'
    write(*,'(A,F10.2,A,F10.2,A)') 'Wet Stage / Base Dry Effic.  = ', eta_wet_actual*100.0d0, ' % / ', eta_dry_target*100.0d0, ' %'
    write(*,'(A,F10.2,A,F10.3,A)') 'Specific Enthalpy Drop (Δh)  = ', delta_h_kJkg, ' kJ/kg (Reheat: ', reheat_factor, ')'
    write(*,'(A,F10.1,A,F10.1,A)') 'Inlet / Exit Enthalpy h1, h2 = ', h1_kJkg, ' kJ/kg / ', h2_kJkg, ' kJ/kg'
    write(*,'(A)') '============================================================'

end program multistage_steam_turbine_wilson_loss