Binary Diffusivity Estimator

program binary_diffusivity
    implicit none
    integer :: mode, i, n_points, iostat_val
    double precision :: T_C, T_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB
    double precision :: VA, mu_cP, phi, rho, D_cm2s, D_m2s, Sc
    double precision :: sigmaAB, epsAB, Tstar, OmegaD, Tprof_C, Tprof_K, Dprof_cm2s, Dprof_m2s
    character(len=80) :: method_name

    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid model selection input.'
        stop
    end if
    read(*,*,iostat=iostat_val) T_C
    read(*,*,iostat=iostat_val) P_atm
    read(*,*,iostat=iostat_val) MA
    read(*,*,iostat=iostat_val) MB
    read(*,*,iostat=iostat_val) sigmaA
    read(*,*,iostat=iostat_val) sigmaB
    read(*,*,iostat=iostat_val) epsA
    read(*,*,iostat=iostat_val) epsB
    read(*,*,iostat=iostat_val) VA
    read(*,*,iostat=iostat_val) mu_cP
    read(*,*,iostat=iostat_val) phi
    read(*,*,iostat=iostat_val) rho
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all binary diffusivity parameters.'
        stop
    end if

    T_K = T_C + 273.15d0
    if (T_K <= 0.0d0) then
        write(*,*) 'ERROR: Temperature must be above absolute zero.'
        stop
    end if
    if (MA <= 0.0d0 .or. MB <= 0.0d0) then
        write(*,*) 'ERROR: Molecular weights must be positive.'
        stop
    end if

    if (mode == 1) then
        if (P_atm <= 0.0d0 .or. sigmaA <= 0.0d0 .or. sigmaB <= 0.0d0 .or. epsA <= 0.0d0 .or. epsB <= 0.0d0) then
            write(*,*) 'ERROR: Gas pressure, sigma, and epsilon/k inputs must be positive.'
            stop
        end if
        method_name = 'Chapman-Enskog gas diffusion'
        call gas_diffusivity(T_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB, D_cm2s, sigmaAB, epsAB, Tstar, OmegaD)
        D_m2s = D_cm2s * 1.0d-4
        Sc = 0.0d0
    else if (mode == 2) then
        if (VA <= 0.0d0 .or. mu_cP <= 0.0d0 .or. phi <= 0.0d0) then
            write(*,*) 'ERROR: Liquid molar volume, viscosity, and association factor must be positive.'
            stop
        end if
        method_name = 'Wilke-Chang liquid diffusion'
        D_cm2s = 7.4d-8 * sqrt(phi * MB) * T_K / (mu_cP * VA**0.6d0)
        D_m2s = D_cm2s * 1.0d-4
        if (rho > 0.0d0) then
            Sc = (mu_cP * 1.0d-3) / (rho * D_m2s)
        else
            Sc = 0.0d0
        end if
        sigmaAB = 0.0d0; epsAB = 0.0d0; Tstar = 0.0d0; OmegaD = 0.0d0
    else
        write(*,*) 'ERROR: Invalid model. Use 1 for gas or 2 for liquid.'
        stop
    end if

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   BINARY DIFFUSIVITY ESTIMATOR ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Method                  = ', trim(method_name)
    write(*,'(A,F12.4,A)')  '  Temperature             = ', T_C, ' deg-C'
    write(*,'(A,F12.4,A)')  '  Temperature (K)         = ', T_K, ' K'
    write(*,'(A,F12.6,A)')  '  Pressure (atm)          = ', P_atm, ' atm'
    write(*,*)
    write(*,'(A)') '--- SPECIES / MOLECULAR INPUTS -----------------------------'
    write(*,'(A,F12.4,A)')  '  Molecular Weight A      = ', MA, ' g/mol'
    write(*,'(A,F12.4,A)')  '  Molecular Weight B      = ', MB, ' g/mol'
    if (mode == 1) then
        write(*,'(A,F12.4,A)') '  sigma_A                 = ', sigmaA, ' Angstrom'
        write(*,'(A,F12.4,A)') '  sigma_B                 = ', sigmaB, ' Angstrom'
        write(*,'(A,F12.4,A)') '  epsilon/k_A             = ', epsA, ' K'
        write(*,'(A,F12.4,A)') '  epsilon/k_B             = ', epsB, ' K'
        write(*,'(A,F12.4,A)') '  sigma_AB                = ', sigmaAB, ' Angstrom'
        write(*,'(A,F12.4,A)') '  epsilon/k_AB            = ', epsAB, ' K'
        write(*,'(A,F12.4)')   '  Reduced Temperature T*  = ', Tstar
        write(*,'(A,F12.4)')   '  Collision Integral      = ', OmegaD
    else
        write(*,'(A,F12.4,A)') '  Solute Molar Volume VA  = ', VA, ' cm3/mol'
        write(*,'(A,F12.4,A)') '  Solvent Viscosity       = ', mu_cP, ' cP'
        write(*,'(A,F12.4)')   '  Association Factor phi  = ', phi
        write(*,'(A,F12.4,A)') '  Liquid Density          = ', rho, ' kg/m3'
    end if
    write(*,*)
    write(*,'(A)') '--- DIFFUSIVITY RESULTS -------------------------------------'
    write(*,'(A,ES12.4,A)') '  Binary Diffusivity (D_AB) = ', D_m2s, ' m2/s'
    write(*,'(A,ES12.4,A)') '  D_AB [cm2/s]             = ', D_cm2s, ' cm2/s'
    if (mode == 2 .and. Sc > 0.0d0) write(*,'(A,ES12.4)') '  Schmidt Number (Sc)      = ', Sc
    write(*,*)
    write(*,'(A)') '--- PARAMETRIC PROFILE --------------------------------------'
    write(*,'(A)') '  T_C [C]      P [atm]       D_AB [m2/s]    D_AB [cm2/s]'
    write(*,'(A)') '  ----------------------------------------------------------'
    n_points = 40
    do i = 1, n_points
        Tprof_C = max(-250.0d0, T_C - 50.0d0 + 100.0d0 * dble(i-1) / dble(n_points-1))
        Tprof_K = Tprof_C + 273.15d0
        if (mode == 1) then
            call gas_diffusivity(Tprof_K, P_atm, MA, MB, sigmaA, sigmaB, epsA, epsB, Dprof_cm2s, sigmaAB, epsAB, Tstar, OmegaD)
        else
            Dprof_cm2s = 7.4d-8 * sqrt(phi * MB) * Tprof_K / (mu_cP * VA**0.6d0)
        end if
        Dprof_m2s = Dprof_cm2s * 1.0d-4
        write(*,'(F9.3,2X,F10.5,2X,ES12.4,2X,ES12.4)') Tprof_C, P_atm, Dprof_m2s, Dprof_cm2s
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    if (mode == 1) then
        write(*,'(A)') '  Chapman-Enskog: D_AB = 0.001858 T^(3/2) sqrt(1/MA+1/MB)/(P sigma_AB^2 Omega_D) [cm2/s]'
        write(*,'(A)') '  sigma_AB = (sigma_A + sigma_B)/2; epsilon_AB/k = sqrt(epsilon_A/k * epsilon_B/k)'
        write(*,'(A)') '  Omega_D uses the Neufeld et al. collision-integral fit.'
    else
        write(*,'(A)') '  Wilke-Chang: D_AB = 7.4E-8 sqrt(phi MB) T/(mu_B VA^0.6) [cm2/s]'
        write(*,'(A)') '  Intended for dilute solutes diffusing in liquid solvents.'
    end if

contains
    subroutine gas_diffusivity(TK, Patm, M1, M2, sig1, sig2, eps1, eps2, Dcm, sigAB, epsAB, Tst, Om)
        implicit none
        double precision, intent(in) :: TK, Patm, M1, M2, sig1, sig2, eps1, eps2
        double precision, intent(out) :: Dcm, sigAB, epsAB, Tst, Om
        sigAB = 0.5d0 * (sig1 + sig2)
        epsAB = sqrt(eps1 * eps2)
        Tst = TK / epsAB
        Om = collision_integral(Tst)
        Dcm = 0.001858d0 * TK**1.5d0 * sqrt(1.0d0/M1 + 1.0d0/M2) / (Patm * sigAB**2 * Om)
    end subroutine gas_diffusivity

    double precision function collision_integral(Ts)
        implicit none
        double precision, intent(in) :: Ts
        collision_integral = 1.06036d0 / Ts**0.15610d0 + 0.19300d0 / exp(0.47635d0*Ts) + &
                             1.03587d0 / exp(1.52996d0*Ts) + 1.76474d0 / exp(3.89411d0*Ts)
    end function collision_integral
end program binary_diffusivity