Gas Permeation Membrane Separation

program gas_permeation_membrane
    implicit none
    integer :: iostat_val, n_elem, i
    double precision :: Q_A_GPU, alpha_AB, Ph_bar, Pl_bar, F_Nm3h, xf_A, stage_cut
    double precision :: Q_A_SI, Q_B_SI, phi_press, xp, yp, xr, Am_total_m2
    double precision :: F_mol_s, P_mol_s, R_mol_s, rec_A_pct, rec_B_pct
    double precision :: dx, x_local, y_local, dAm, dP_Pa, Ph_Pa, Pl_Pa
    double precision, parameter :: GPU_TO_SI = 3.348d-10 ! mol / (m2.s.Pa)

    ! Read inputs
    read(*,*,iostat=iostat_val) Q_A_GPU       ! Fast Gas Permeance [GPU] (e.g. 100.0)
    read(*,*,iostat=iostat_val) alpha_AB      ! Selectivity alpha (QA / QB) (e.g. 25.0)
    read(*,*,iostat=iostat_val) Ph_bar        ! Feed Shell High Pressure [bar a] (e.g. 30.0)
    read(*,*,iostat=iostat_val) Pl_bar        ! Permeate Lumen Low Pressure [bar a] (e.g. 1.2)
    read(*,*,iostat=iostat_val) F_Nm3h        ! Feed Flow Rate [Nm3/h] (e.g. 500.0)
    read(*,*,iostat=iostat_val) xf_A          ! Feed Fast Gas Mole Fraction (e.g. 0.40)
    read(*,*,iostat=iostat_val) stage_cut     ! Target Stage Cut theta (Perm / Feed) (e.g. 0.35)

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

    if (alpha_AB <= 1.0d0 .or. Ph_bar <= Pl_bar .or. stage_cut <= 0.0d0 .or. stage_cut >= 1.0d0) then
        write(*,*) 'ERROR: alpha must be > 1, Ph > Pl, and 0 < stage_cut < 1.'
        stop
    end if

    Ph_Pa = Ph_bar * 1.0d5
    Pl_Pa = Pl_bar * 1.0d5
    phi_press = Pl_bar / Ph_bar
    Q_A_SI = Q_A_GPU * GPU_TO_SI
    Q_B_SI = (Q_A_GPU / alpha_AB) * GPU_TO_SI

    ! Molar flows (1 Nm3 = 44.615 mol)
    F_mol_s = (F_Nm3h / 3600.0d0) * 44.615d0
    P_mol_s = stage_cut * F_mol_s
    R_mol_s = F_mol_s - P_mol_s

    ! Discretized cross-flow permeation model across membrane length
    n_elem = 200
    Am_total_m2 = 0.0d0
    x_local = xf_A
    yp = 0.0d0

    ! Solve average permeate concentration yp
    ! For binary cross flow Weller-Steiner / Shindo approximate solution
    yp = (1.0d0 / (2.0d0 * (1.0d0 - 1.0d0/alpha_AB))) * &
         (1.0d0 + (alpha_AB - 1.0d0)*phi_press + (alpha_AB - 1.0d0)*xf_A - &
         sqrt((1.0d0 + (alpha_AB - 1.0d0)*phi_press + (alpha_AB - 1.0d0)*xf_A)**2 - &
         4.0d0 * alpha_AB * (alpha_AB - 1.0d0) * phi_press * xf_A))
    if (yp <= 0.0d0 .or. yp >= 1.0d0) yp = xf_A * 1.8d0
    if (yp > 0.999d0) yp = 0.999d0

    ! Retentate concentration xr by overall component balance
    xr = (F_mol_s * xf_A - P_mol_s * yp) / R_mol_s
    if (xr < 0.001d0) xr = 0.001d0

    ! Membrane Area required
    ! Log-mean partial pressure driving force
    Am_total_m2 = (P_mol_s * yp) / max(1.0d-4, (Q_A_SI * (Ph_Pa * 0.5d0 * (xf_A + xr) - Pl_Pa * yp)))

    rec_A_pct = (P_mol_s * yp) / (F_mol_s * xf_A) * 100.0d0
    rec_B_pct = (R_mol_s * (1.0d0 - xr)) / (F_mol_s * (1.0d0 - xf_A)) * 100.0d0

    ! Output Results
    write(*,'(A)') '============================================================'
    write(*,'(A)') ' THERMOFLUIDCALC — GAS PERMEATION MEMBRANE SEPARATION'
    write(*,'(A)') '============================================================'
    write(*,'(A,F10.1,A,F10.2)') 'Fast Gas Permeance / Selectivity  = ', Q_A_GPU, ' GPU / ', alpha_AB
    write(*,'(A,F10.1,A,F10.1,A)') 'Feed / Permeate Pressures         = ', Ph_bar, ' bar / ', Pl_bar, ' bar'
    write(*,'(A,F10.1,A,F6.3)')  'Feed Flow Rate / Mole Fraction    = ', F_Nm3h, ' Nm3/h / ', xf_A
    write(*,'(A)') '------------------------------------------------------------'
    write(*,'(A,F10.2,A)')  'Target Stage Cut (theta)          = ', stage_cut * 100.0d0, ' %'
    write(*,'(A,F10.3,A,F10.3)') 'Permeate Purity (yp) / Retentate  = ', yp, ' / ', xr
    write(*,'(A,F10.1,A)')  'Required Membrane Surface Area    = ', Am_total_m2, ' m2'
    write(*,'(A,F10.1,A,F10.1,A)') 'Permeate / Retentate Flow Rates   = ', &
        P_mol_s * 3600.0d0 / 44.615d0, ' Nm3/h / ', R_mol_s * 3600.0d0 / 44.615d0, ' Nm3/h'
    write(*,'(A,F10.1,A,F10.1,A)') 'Fast Gas Recovery / Slow Gas Rec  = ', rec_A_pct, ' % / ', rec_B_pct, ' %'
    write(*,'(A)') '============================================================'

end program gas_permeation_membrane