Counter-Current Extraction (Kremser)

program countercurrent_extraction_kremser
    implicit none
    integer :: iostat_val
    double precision :: F_kgh, S_kgh, xF_pct, yS_pct, target_rec_pct, m_KD, HTU_m
    double precision :: xF, yS, xN, y1, E_factor, N_theor, NTU_OL, Z_pack_m
    double precision :: Extracted_kgh, Solute_in_kgh, Solute_out_kgh

    ! Read inputs
    read(*,*,iostat=iostat_val) F_kgh           ! Feed Flow Rate [kg/h] (e.g. 1000.0)
    read(*,*,iostat=iostat_val) S_kgh           ! Solvent Flow Rate [kg/h] (e.g. 800.0)
    read(*,*,iostat=iostat_val) xF_pct          ! Feed Solute Weight [%] (e.g. 10.0)
    read(*,*,iostat=iostat_val) yS_pct          ! Solvent Solute Content [%] (e.g. 0.0)
    read(*,*,iostat=iostat_val) target_rec_pct  ! Target Recovery [%] (e.g. 95.0)
    read(*,*,iostat=iostat_val) m_KD            ! Distribution Coefficient KD = y/x (e.g. 1.80)
    read(*,*,iostat=iostat_val) HTU_m           ! Height of Transfer Unit HTU [m] (e.g. 0.75)

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

    if (F_kgh <= 0.0d0 .or. S_kgh <= 0.0d0 .or. m_KD <= 0.0d0 .or. target_rec_pct >= 100.0d0) then
        write(*,*) 'ERROR: Flows and KD must be positive, recovery < 100%.'
        stop
    end if

    xF = xF_pct / 100.0d0
    yS = yS_pct / 100.0d0

    ! Target raffinate concentration
    xN = xF * (1.0d0 - target_rec_pct / 100.0d0)

    ! Extraction Factor E = m * S / F
    E_factor = (m_KD * S_kgh) / F_kgh

    ! Kremser equation for N_theor
    if (abs(E_factor - 1.0d0) > 0.01d0) then
        N_theor = log(((xF - yS/m_KD) / (xN - yS/m_KD)) * (1.0d0 - 1.0d0/E_factor) + 1.0d0/E_factor) / log(E_factor)
        NTU_OL = (E_factor / (E_factor - 1.0d0)) * &
                 log(((1.0d0 - 1.0d0/E_factor) * ((xF - yS/m_KD) / (xN - yS/m_KD))) + 1.0d0/E_factor)
    else
        N_theor = (xF - xN) / max(1.0d-5, (xN - yS/m_KD))
        NTU_OL = N_theor
    end if

    if (N_theor < 0.1d0) N_theor = 0.5d0
    if (NTU_OL < 0.1d0) NTU_OL = 0.5d0

    Z_pack_m = NTU_OL * HTU_m

    Solute_in_kgh = F_kgh * xF
    Extracted_kgh = Solute_in_kgh * (target_rec_pct / 100.0d0)
    Solute_out_kgh = Solute_in_kgh - Extracted_kgh
    y1 = (Extracted_kgh + S_kgh * yS) / S_kgh

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — COUNTER-CURRENT EXTRACTION (KREMSER)'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.1,A)') 'Feed / Solvent Flow Rate  = ', F_kgh, ' kg/h / ', S_kgh, ' kg/h'
    write(*,'(A,F10.2,A,F10.2)') 'Feed Solute / KD (y/x)    = ', xF_pct, ' % / ', m_KD
    write(*,'(A,F10.3)')    'Extraction Factor (E)     = ', E_factor
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Theoretical Stages (N)    = ', N_theor, ' stages'
    write(*,'(A,F10.2,A)')  'Number of Transfer Units  = ', NTU_OL, ' NTU'
    write(*,'(A,F10.2,A)')  'Packed Column Height (Z)  = ', Z_pack_m, ' meters (HTU = ', HTU_m, ' m)'
    write(*,'(A,F10.2,A)')  'Solute Recovery Yield     = ', target_rec_pct, ' %'
    write(*,'(A,F10.3,A,F10.3,A)') 'Final Raffinate / Extract = ', xN*100.0d0, ' % / ', y1*100.0d0, ' %'
    write(*,'(A,F10.2,A)')  'Mass Solute Extracted     = ', Extracted_kgh, ' kg/h'
    write(*,'(A)') '============================================================'

end program countercurrent_extraction_kremser