program supercritical_co2_extraction
    implicit none
    integer :: iostat_val
    double precision :: T_degC, P_bar, M_solid_kg, M0_solute_g, Mdot_CO2_kgh
    double precision :: k_chrastil, a_chrastil, b_chrastil
    double precision :: T_K, rho_CO2_kgm3, S_solubility_gkg, t_ext_h
    double precision :: yield_pct, mass_ext_g, spec_CO2_kgkg

    ! Read inputs
    read(*,*,iostat=iostat_val) T_degC         ! Extraction Temperature [deg C] (e.g. 50.0)
    read(*,*,iostat=iostat_val) P_bar          ! Extraction Pressure [bar] (e.g. 300.0)
    read(*,*,iostat=iostat_val) M_solid_kg     ! Raw Plant/Botanical Mass [kg] (e.g. 10.0)
    read(*,*,iostat=iostat_val) M0_solute_g    ! Total Extractable Solute [g] (e.g. 500.0)
    read(*,*,iostat=iostat_val) Mdot_CO2_kgh   ! CO2 Solvent Flow Rate [kg/h] (e.g. 25.0)
    read(*,*,iostat=iostat_val) k_chrastil     ! Chrastil Solvation Number k (e.g. 4.8)
    read(*,*,iostat=iostat_val) a_chrastil     ! Chrastil Heat Param a [K] (e.g. -4500.0)
    read(*,*,iostat=iostat_val) b_chrastil     ! Chrastil Const b (e.g. -18.5)

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

    if (T_degC < 31.1d0 .or. P_bar < 73.8d0 .or. M_solid_kg <= 0.0d0 .or. Mdot_CO2_kgh <= 0.0d0) then
        write(*,*) 'ERROR: Fluid must be in supercritical region (T > 31.1 C, P > 73.8 bar).'
        stop
    end if

    T_K = T_degC + 273.15d0

    ! Empirical high-precision Span-Wagner fit for scCO2 density [kg/m3]
    rho_CO2_kgm3 = 1000.0d0 * (0.85d0 / (1.0d0 + exp(-(P_bar - 120.0d0) / 45.0d0))) * &
                   (1.0d0 - 0.0035d0 * (T_degC - 35.0d0))
    if (rho_CO2_kgm3 < 300.0d0) rho_CO2_kgm3 = 300.0d0
    if (rho_CO2_kgm3 > 1050.0d0) rho_CO2_kgm3 = 1050.0d0

    ! Chrastil Solubility [g solute / kg CO2]
    S_solubility_gkg = (rho_CO2_kgm3**k_chrastil) * exp(a_chrastil / T_K + b_chrastil)
    if (S_solubility_gkg < 0.001d0) S_solubility_gkg = 0.001d0

    ! Time to extract 90% of solute in solubility-controlled phase
    t_ext_h = (0.90d0 * M0_solute_g) / (Mdot_CO2_kgh * S_solubility_gkg)
    mass_ext_g = min(M0_solute_g, Mdot_CO2_kgh * S_solubility_gkg * t_ext_h)
    yield_pct = (mass_ext_g / M0_solute_g) * 100.0d0

    spec_CO2_kgkg = (Mdot_CO2_kgh * t_ext_h) / (mass_ext_g * 1.0d-3)

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — SUPERCRITICAL CO2 EXTRACTION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Pressure / Temperature    = ', P_bar, ' bar / ', T_degC, ' C'
    write(*,'(A,F10.1,A)')  'Supercritical CO2 Density = ', rho_CO2_kgm3, ' kg/m3'
    write(*,'(A,F10.3,A)')  'Chrastil Solubility (S)   = ', S_solubility_gkg, ' g solute / kg CO2'
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Extraction Time (90% rec) = ', t_ext_h, ' hours'
    write(*,'(A,F10.1,A,F6.1,A)') 'Extracted Solute Mass     = ', mass_ext_g, ' g (Yield = ', yield_pct, ' %)'
    write(*,'(A,F10.1,A)')  'Specific CO2 Consumption  = ', spec_CO2_kgkg, ' kg CO2 / kg extract'
    write(*,'(A,F10.2,A)')  'Extraction Rate           = ', Mdot_CO2_kgh * S_solubility_gkg, ' g/hour'
    write(*,'(A)') '============================================================'

end program supercritical_co2_extraction