program mass_diffusion
    implicit none

    ! Inputs
    integer :: mode ! 1 = Equimolar Counter-Diffusion (ECD), 2 = Diffusion in Stagnant Gas (Stefan)
    double precision :: T_C, T_K ! Temperature in C and Kelvin
    double precision :: P_kPa, P_Pa ! Pressure in kPa and Pascals
    double precision :: D_AB ! Diffusivity in m2/s
    double precision :: L ! Film thickness in m
    double precision :: P_A1_kPa, P_A1_Pa ! Boundary 1 partial pressure of solute A in kPa and Pa
    double precision :: P_A2_kPa, P_A2_Pa ! Boundary 2 partial pressure of solute A in kPa and Pa
    double precision :: M_A ! Molecular weight of solute A in g/mol

    ! Constants
    double precision, parameter :: R = 8.314462618d0 ! Universal gas constant, J/(mol.K)

    ! Physical properties
    double precision :: C_total ! Total concentration, mol/m3
    double precision :: y_A1, y_A2 ! Mole fractions at boundaries
    double precision :: P_B1, P_B2, P_Blm ! Partial pressures of stagnant gas B
    double precision :: N_A ! Molar flux of A, mol/m2.s
    double precision :: n_A_mass ! Mass flux of A, g/m2.s
    double precision :: h_m ! Mass transfer coefficient, m/s

    ! Local profile variables
    integer :: i, n_points
    double precision :: x, dx, local_yA, local_CA, local_PA
    character(len=50) :: mode_name
    integer :: iostat_val

    ! Read inputs
    read(*,*,iostat=iostat_val) mode
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid mode input (must be 1 or 2).'
        stop
    end if

    read(*,*,iostat=iostat_val) T_C
    read(*,*,iostat=iostat_val) P_kPa
    read(*,*,iostat=iostat_val) D_AB
    read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) P_A1_kPa
    read(*,*,iostat=iostat_val) P_A2_kPa
    read(*,*,iostat=iostat_val) M_A

    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all mass diffusion parameters.'
        stop
    end if

    ! Validations
    if (mode /= 1 .and. mode /= 2) then
        write(*,*) 'ERROR: Mode must be 1 (Equimolar Counter-Diffusion) or 2 (Stagnant Film).'
        stop
    end if
    if (T_C < -273.15d0) then
        write(*,*) 'ERROR: Temperature must be above absolute zero.'
        stop
    end if
    if (P_kPa <= 0.0d0) then
        write(*,*) 'ERROR: Pressure must be positive.'
        stop
    end if
    if (D_AB <= 0.0d0) then
        write(*,*) 'ERROR: Diffusivity must be positive.'
        stop
    end if
    if (L <= 0.0d0) then
        write(*,*) 'ERROR: Film thickness must be positive.'
        stop
    end if
    if (P_A1_kPa < 0.0d0 .or. P_A1_kPa > P_kPa) then
        write(*,*) 'ERROR: Boundary 1 partial pressure must be between 0 and total pressure.'
        stop
    end if
    if (P_A2_kPa < 0.0d0 .or. P_A2_kPa > P_kPa) then
        write(*,*) 'ERROR: Boundary 2 partial pressure must be between 0 and total pressure.'
        stop
    end if
    if (M_A <= 0.0d0) then
        write(*,*) 'ERROR: Molecular weight must be positive.'
        stop
    end if

    ! Conversions
    T_K = T_C + 273.15d0
    P_Pa = P_kPa * 1000.0d0
    P_A1_Pa = P_A1_kPa * 1000.0d0
    P_A2_Pa = P_A2_kPa * 1000.0d0

    ! Total molar concentration (ideal gas law C = P / RT)
    C_total = P_Pa / (R * T_K)
    y_A1 = P_A1_Pa / P_Pa
    y_A2 = P_A2_Pa / P_Pa

    ! Core Calculations
    if (mode == 1) then
        mode_name = "Equimolar Counter-Diffusion"
        
        ! Molar flux N_A = D_AB * (C_A1 - C_A2) / L
        N_A = D_AB * C_total * (y_A1 - y_A2) / L
        h_m = D_AB / L
        P_Blm = 0.0d0 ! Not applicable
    else
        mode_name = "Diffusion through a Stagnant Gas B"
        
        ! Boundary partial pressures of stagnant B
        P_B1 = P_Pa - P_A1_Pa
        P_B2 = P_Pa - P_A2_Pa
        
        ! Log-mean partial pressure of B
        if (abs(P_B1 - P_B2) < 1.0d-5 * P_Pa) then
            P_Blm = (P_B1 + P_B2) / 2.0d0
        else
            P_Blm = (P_B2 - P_B1) / log(P_B2 / P_B1)
        end if
        
        ! Molar flux N_A = D_AB * P * (P_A1 - P_A2) / (R * T * L * P_Blm)
        N_A = D_AB * P_Pa * (P_A1_Pa - P_A2_Pa) / (R * T_K * L * P_Blm)
        h_m = D_AB * P_Pa / (L * P_Blm)
    end if

    ! Convert to mass flux (g/m2.s)
    n_A_mass = N_A * M_A

    ! Output results
    write(*,'(A)') '============================================================'
    write(*,'(A)') '      MASS DIFFUSION CALCULATION ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Calculation Mode         = ', trim(mode_name)
    write(*,'(A,I2)')         '  Mode Code                = ', mode
    write(*,'(A,F12.2,A)')  '  Temperature              = ', T_C, ' deg-C'
    write(*,'(A,F12.2,A)')  '  Total Pressure           = ', P_kPa, ' kPa'
    write(*,'(A,ES12.4,A)') '  Binary Diffusivity D_AB  = ', D_AB, ' m2/s'
    write(*,'(A,F12.4,A)')  '  Film Thickness (L)       = ', L, ' m'
    write(*,'(A,F12.2,A)')  '  Molecular Weight M_A     = ', M_A, ' g/mol'
    write(*,*)
    write(*,'(A)') '--- BOUNDARY CONDITIONS -------------------------------------'
    write(*,'(A,F12.4,A)')  '  Boundary 1 Partial P_A1  = ', P_A1_kPa, ' kPa'
    write(*,'(A,F12.4,A)')  '  Boundary 2 Partial P_A2  = ', P_A2_kPa, ' kPa'
    write(*,'(A,F12.6)')    '  Boundary 1 Mole Frac yA1 = ', y_A1
    write(*,'(A,F12.6)')    '  Boundary 2 Mole Frac yA2 = ', y_A2
    write(*,*)
    write(*,'(A)') '--- ENGINE RESULTS ------------------------------------------'
    write(*,'(A,F12.4,A)')  '  Total Concentration C    = ', C_total, ' mol/m3'
    if (mode == 2) then
        write(*,'(A,F12.4,A)')  '  Log-Mean Pressure P_Blm  = ', P_Blm/1000.0d0, ' kPa'
    end if
    write(*,'(A,F12.6,A)')  '  Mass Transfer Coeff (hm) = ', h_m, ' m/s'
    write(*,'(A,ES12.4,A)') '  Molar Flux of A (N_A)    = ', N_A, ' mol/m2.s'
    write(*,'(A,ES12.4,A)') '  Mass Flux of A (n_A)     = ', n_A_mass, ' g/m2.s'
    write(*,*)

    ! ============================================
    ! LOCAL PROFILE DATA (along path coordinate x)
    ! ============================================
    write(*,'(A)') '--- LOCAL PROFILE ALONG FILM --------------------------------'
    write(*,'(A)') '  x [m]       y_A [-]       C_A [mol/m3]  P_A [kPa]'
    write(*,'(A)') '  -----------------------------------------------------------'

    n_points = 40
    dx = L / dble(n_points)

    ! Include x = 0 initially
    write(*,'(F8.4,4X,F10.6,4X,F12.4,4X,F10.4)') 0.0d0, y_A1, y_A1 * C_total, P_A1_kPa

    do i = 1, n_points
        x = dble(i) * dx
        
        if (mode == 1) then
            ! Linear profile for ECD
            local_yA = y_A1 - (y_A1 - y_A2) * (x / L)
        else
            ! Non-linear exponential-like profile for stagnant gas
            if (abs(y_A1 - y_A2) < 1.0d-7) then
                local_yA = y_A1
            else
                local_yA = 1.0d0 - (1.0d0 - y_A1) * ((1.0d0 - y_A2)/(1.0d0 - y_A1))**(x/L)
            end if
        end if
        
        local_CA = local_yA * C_total
        local_PA = local_yA * P_kPa
        
        write(*,'(F8.4,4X,F10.6,4X,F12.4,4X,F10.4)') x, local_yA, local_CA, local_PA
    end do

    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    if (mode == 1) then
        write(*,'(A)') "  Fick's First Law for Equimolar Counter-Diffusion (linear gradient)."
    else
        write(*,'(A)') "  Stefan's Diffusion Law (diffusion through stagnant gas B, logarithmic profile)."
    end if
    write(*,'(A)') '  Total concentration computed from ideal gas equation of state.'

end program mass_diffusion
