program membrane_diffusion
    implicit none
    integer :: model, i, n_points, iostat_val
    double precision :: thickness, area, T_C, pFeed_bar, pPerm_bar, yAf, yAp, yBf, yBp
    double precision :: DA, SA, DB, SB, rpore, porosity, tau, MWA, MWB, feedFlow
    double precision :: T, pF, pP, R, permA, permB, permeanceA, permeanceB, JA, JB, JmassA
    double precision :: dPA, dPB, DKA, DKB, DKeffA, DKeffB, PsolA, PsolB, PeffA, PeffB
    double precision :: selectivity, flowA, flowB, stageCut, purityA
    double precision :: pSweep, dPAs, dPBs, JAs, JBs, flowAs, flowBs, thetaS, purityS
    character(len=80) :: model_name

    read(*,*,iostat=iostat_val) model
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid model input.'
        stop
    end if
    read(*,*,iostat=iostat_val) thickness
    read(*,*,iostat=iostat_val) area
    read(*,*,iostat=iostat_val) T_C
    read(*,*,iostat=iostat_val) pFeed_bar
    read(*,*,iostat=iostat_val) pPerm_bar
    read(*,*,iostat=iostat_val) yAf
    read(*,*,iostat=iostat_val) yAp
    read(*,*,iostat=iostat_val) yBf
    read(*,*,iostat=iostat_val) yBp
    read(*,*,iostat=iostat_val) DA
    read(*,*,iostat=iostat_val) SA
    read(*,*,iostat=iostat_val) DB
    read(*,*,iostat=iostat_val) SB
    read(*,*,iostat=iostat_val) rpore
    read(*,*,iostat=iostat_val) porosity
    read(*,*,iostat=iostat_val) tau
    read(*,*,iostat=iostat_val) MWA
    read(*,*,iostat=iostat_val) MWB
    read(*,*,iostat=iostat_val) feedFlow
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all membrane diffusion inputs.'
        stop
    end if

    if (thickness <= 0.0d0 .or. area <= 0.0d0 .or. pFeed_bar < 0.0d0 .or. pPerm_bar < 0.0d0) then
        write(*,*) 'ERROR: Thickness, area, and pressures must be valid positive values.'
        stop
    end if
    if (MWA <= 0.0d0 .or. MWB <= 0.0d0) then
        write(*,*) 'ERROR: Molecular weights must be positive.'
        stop
    end if
    if (T_C <= -273.15d0) then
        write(*,*) 'ERROR: Temperature must be above absolute zero.'
        stop
    end if

    R = 8.314462618d0
    T = T_C + 273.15d0
    pF = pFeed_bar * 1.0d5
    pP = pPerm_bar * 1.0d5
    dPA = yAf*pF - yAp*pP
    dPB = yBf*pF - yBp*pP

    PsolA = DA*SA
    PsolB = DB*SB
    DKA = (2.0d0/3.0d0)*rpore*sqrt(8.0d0*R*T/(3.141592653589793d0*(MWA/1000.0d0)))
    DKB = (2.0d0/3.0d0)*rpore*sqrt(8.0d0*R*T/(3.141592653589793d0*(MWB/1000.0d0)))
    if (tau > 0.0d0) then
        DKeffA = porosity/tau*DKA
        DKeffB = porosity/tau*DKB
    else
        DKeffA = 0.0d0
        DKeffB = 0.0d0
    end if

    select case(model)
    case(1)
        model_name = 'Solution-diffusion dense membrane'
        PeffA = PsolA
        PeffB = PsolB
        permeanceA = PeffA/thickness
        permeanceB = PeffB/thickness
        JA = permeanceA*dPA
        JB = permeanceB*dPB
    case(2)
        model_name = 'Knudsen porous membrane'
        PeffA = DKeffA/(R*T)
        PeffB = DKeffB/(R*T)
        permeanceA = PeffA/thickness
        permeanceB = PeffB/thickness
        JA = permeanceA*dPA
        JB = permeanceB*dPB
    case(3)
        model_name = 'Mixed solution-diffusion plus Knudsen resistance'
        PeffA = series_perm(PsolA, DKeffA/(R*T))
        PeffB = series_perm(PsolB, DKeffB/(R*T))
        permeanceA = PeffA/thickness
        permeanceB = PeffB/thickness
        JA = permeanceA*dPA
        JB = permeanceB*dPB
    case default
        write(*,*) 'ERROR: Invalid model. Use 1 solution, 2 Knudsen, 3 mixed.'
        stop
    end select

    if (permeanceB > 0.0d0) then
        selectivity = permeanceA/permeanceB
    else
        selectivity = 0.0d0
    end if
    flowA = JA*area
    flowB = JB*area
    JmassA = JA*MWA/1000.0d0
    if (feedFlow > 0.0d0) then
        stageCut = (flowA + flowB)/feedFlow
    else
        stageCut = 0.0d0
    end if
    if ((flowA + flowB) > 0.0d0) then
        purityA = flowA/(flowA + flowB)
    else
        purityA = 0.0d0
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   MEMBRANE DIFFUSION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Model                    = ', trim(model_name)
    write(*,'(A,ES12.4,A)') '  Thickness                = ', thickness, ' m'
    write(*,'(A,ES12.4,A)') '  Membrane Area            = ', area, ' m2'
    write(*,'(A,F12.4,A)')  '  Temperature              = ', T_C, ' deg-C'
    write(*,'(A,ES12.4,A)') '  Feed Pressure            = ', pFeed_bar, ' bar'
    write(*,'(A,ES12.4,A)') '  Permeate Pressure        = ', pPerm_bar, ' bar'
    write(*,*)
    write(*,'(A)') '--- TRANSPORT PROPERTIES ------------------------------------'
    write(*,'(A,ES12.4,A)') '  Solution D A             = ', DA, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Solution D B             = ', DB, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Solubility A             = ', SA, ' mol/m3.Pa'
    write(*,'(A,ES12.4,A)') '  Solubility B             = ', SB, ' mol/m3.Pa'
    write(*,'(A,ES12.4,A)') '  Knudsen D A              = ', DKA, ' m2/s'
    write(*,'(A,ES12.4,A)') '  Knudsen D B              = ', DKB, ' m2/s'
    write(*,'(A,ES12.4)')   '  Porosity/Tortuosity      = ', merge(porosity/tau, 0.0d0, tau>0.0d0)
    write(*,*)
    write(*,'(A)') '--- MEMBRANE RESULTS ----------------------------------------'
    write(*,'(A,ES12.4,A)') '  Effective Permeability A = ', PeffA, ' mol.m/m2.s.Pa'
    write(*,'(A,ES12.4,A)') '  Effective Permeability B = ', PeffB, ' mol.m/m2.s.Pa'
    write(*,'(A,ES12.4,A)') '  Permeance A              = ', permeanceA, ' mol/m2.s.Pa'
    write(*,'(A,ES12.4,A)') '  Permeance B              = ', permeanceB, ' mol/m2.s.Pa'
    write(*,'(A,ES12.4)')   '  Ideal Selectivity        = ', selectivity
    write(*,'(A,ES12.4,A)') '  Flux A                   = ', JA, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Flux B                   = ', JB, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Mass Flux A              = ', JmassA, ' kg/m2.s'
    write(*,'(A,ES12.4,A)') '  Permeate Flow A          = ', flowA, ' mol/s'
    write(*,'(A,ES12.4,A)') '  Permeate Flow B          = ', flowB, ' mol/s'
    write(*,'(A,ES12.4)')   '  Stage Cut Estimate       = ', stageCut
    write(*,'(A,ES12.4)')   '  Permeate Purity A        = ', purityA
    write(*,*)

    write(*,'(A)') '--- PRESSURE SWEEP PROFILE ----------------------------------'
    write(*,'(A)') '  pFeed[bar]   J_A[mol/m2.s] J_B[mol/m2.s] selectivity   stage_cut    purity_A'
    write(*,'(A)') '  ---------------------------------------------------------------------------'
    n_points = 45
    do i=0,n_points-1
        pSweep = max(pPerm_bar + 0.05d0, 0.5d0*pFeed_bar + 1.5d0*pFeed_bar*dble(i)/dble(n_points-1))
        dPAs = yAf*pSweep*1.0d5 - yAp*pP
        dPBs = yBf*pSweep*1.0d5 - yBp*pP
        JAs = permeanceA*dPAs
        JBs = permeanceB*dPBs
        flowAs = JAs*area
        flowBs = JBs*area
        if (feedFlow > 0.0d0) then
            thetaS = (flowAs + flowBs)/feedFlow
        else
            thetaS = 0.0d0
        end if
        if ((flowAs + flowBs) > 0.0d0) then
            purityS = flowAs/(flowAs + flowBs)
        else
            purityS = 0.0d0
        end if
        write(*,'(F10.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
            pSweep, JAs, JBs, selectivity, thetaS, purityS
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Solution-diffusion: P = D*S; permeance = P/thickness; J = permeance*Delta p.'
    write(*,'(A)') '  Knudsen: D_K = (2/3) r_p sqrt(8RT/(pi MW)); effective D_K = eps/tau D_K.'
    write(*,'(A)') '  Mixed model combines solution and Knudsen permeabilities as series resistances.'

contains
    double precision function series_perm(P1, P2)
        implicit none
        double precision, intent(in) :: P1, P2
        if (P1 <= 0.0d0) then
            series_perm = P2
        else if (P2 <= 0.0d0) then
            series_perm = P1
        else
            series_perm = 1.0d0/(1.0d0/P1 + 1.0d0/P2)
        end if
    end function series_perm
end program membrane_diffusion
