program supercritical_organic_rankine_sorc
    implicit none
    integer :: iostat_val, fluid_type
    double precision :: P_high_bar, T_tit_C, T_cond_C, m_dot_fluid_kgs, eta_t, eta_p, eps_recup
    double precision :: P_crit_bar, T_crit_C, cp_liq, cp_vap, rho_liq, hfg_kJkg
    double precision :: h1, h2, h2s, h3, h4, h5, h6, w_turb, w_pump, w_net, q_in, eta_th, P_net_kW
    double precision :: T_exp_out_C, q_recup_kW

    ! Read inputs
    read(*,*,iostat=iostat_val) fluid_type      ! 1=R245fa, 2=R1233zd(E), 3=n-Pentane, 4=Toluene
    read(*,*,iostat=iostat_val) P_high_bar      ! Supercritical Turbine Inlet Pressure [bar] (e.g. 45.0)
    read(*,*,iostat=iostat_val) T_tit_C         ! Turbine Inlet Temp [deg C] (e.g. 190.0)
    read(*,*,iostat=iostat_val) T_cond_C        ! Condensing Temp [deg C] (e.g. 35.0)
    read(*,*,iostat=iostat_val) m_dot_fluid_kgs ! Working Fluid Mass Flow [kg/s] (e.g. 15.0)
    read(*,*,iostat=iostat_val) eta_t           ! Turbine Isentropic Efficiency (e.g. 0.85)
    read(*,*,iostat=iostat_val) eta_p           ! Pump Isentropic Efficiency (e.g. 0.75)
    read(*,*,iostat=iostat_val) eps_recup       ! Recuperator Effectiveness (e.g. 0.75)

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

    ! Fluid thermophysical properties
    if (fluid_type == 1) then ! R245fa
        P_crit_bar = 36.51d0; T_crit_C = 154.0d0; cp_liq = 1.35d0; cp_vap = 0.95d0; rho_liq = 1320.0d0; hfg_kJkg = 195.0d0
    else if (fluid_type == 2) then ! R1233zd(E)
        P_crit_bar = 35.70d0; T_crit_C = 165.6d0; cp_liq = 1.25d0; cp_vap = 0.90d0; rho_liq = 1260.0d0; hfg_kJkg = 190.0d0
    else if (fluid_type == 3) then ! n-Pentane
        P_crit_bar = 33.70d0; T_crit_C = 196.5d0; cp_liq = 2.40d0; cp_vap = 1.80d0; rho_liq = 610.0d0;  hfg_kJkg = 355.0d0
    else ! Toluene
        P_crit_bar = 41.26d0; T_crit_C = 318.6d0; cp_liq = 1.85d0; cp_vap = 1.45d0; rho_liq = 840.0d0;  hfg_kJkg = 360.0d0
    end if

    if (P_high_bar <= P_crit_bar) then
        write(*,*) 'WARNING: P_high must be supercritical (P > P_crit). Adjusting to 1.15 * P_crit.'
        P_high_bar = 1.15d0 * P_crit_bar
    end if

    ! State 1: Turbine Inlet (Supercritical Vapor)
    h1 = hfg_kJkg + cp_vap * (T_tit_C - T_crit_C) + 400.0d0

    ! State 2: Expander Outlet (Superheated Dry Vapor)
    w_turb = eta_t * (cp_vap * (T_tit_C - T_cond_C) + 0.08d0 * (P_high_bar - 2.5d0))
    h2 = h1 - w_turb
    T_exp_out_C = T_cond_C + (h2 - (hfg_kJkg + 400.0d0 - 50.0d0)) / cp_vap
    if (T_exp_out_C < T_cond_C + 5.0d0) T_exp_out_C = T_cond_C + 5.0d0

    ! State 3: Liquid Condenser Outlet
    h3 = 100.0d0

    ! State 4: High Pressure Pump Outlet
    w_pump = ((P_high_bar - 2.5d0) * 100.0d0 / (rho_liq * eta_p))
    h4 = h3 + w_pump

    ! Recuperation: State 5 (Preheated fluid) & State 6 (Cooled exhaust vapor)
    q_recup_kW = eps_recup * m_dot_fluid_kgs * cp_vap * (T_exp_out_C - (T_cond_C + 2.0d0))
    h5 = h4 + (q_recup_kW / m_dot_fluid_kgs)

    ! Supercritical Heater Heat Input
    q_in = h1 - h5
    w_net = w_turb - w_pump
    eta_th = (w_net / q_in) * 100.0d0
    P_net_kW = m_dot_fluid_kgs * w_net

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERCRITICAL ORGANIC RANKINE CYCLE (sORC)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.2,A,F10.2,A)') 'P_high / P_crit           = ', P_high_bar, ' bar / ', P_crit_bar, ' bar'
    write(*,'(A,F10.1,A,F10.1,A)') 'TIT Temp / Condensing T   = ', T_tit_C, ' deg C / ', T_cond_C, ' deg C'
    write(*,'(A,F10.2,A,F10.2)')   'Mass Flow / Recup Effic.  = ', m_dot_fluid_kgs, ' kg/s / ', eps_recup
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'NET POWER OUTPUT (W_net)  = ', P_net_kW, ' kW'
    write(*,'(A,F10.2,A)')  'THERMAL EFFICIENCY (eta)  = ', eta_th, ' %'
    write(*,'(A,F10.2,A,F10.2,A)') 'Turbine / Pump Power      = ', m_dot_fluid_kgs*w_turb, ' kW / ', m_dot_fluid_kgs*w_pump, ' kW'
    write(*,'(A,F10.2,A)')  'Heater Heat Duty (Q_in)   = ', m_dot_fluid_kgs*q_in, ' kW'
    write(*,'(A,F10.2,A)')  'Recuperator Heat Duty     = ', q_recup_kW, ' kW'
    write(*,'(A,F10.1,A)')  'Turbine Discharge Temp    = ', T_exp_out_C, ' deg C'
    write(*,'(A)') '============================================================'

end program supercritical_organic_rankine_sorc
