      SUBROUTINE INITAL
C...
C...  Three Component Diffusion
C...
C...  Consided the following system consisting of two bulbs connected
C...  by a tube
C...
C...
C...          ........                                  ........
C...          .      .+---------------ZL---------------+.      .
C...          .      ....................................      .
C...          .  VO  ....................................  VL  .
C...          .      .               tube               .      .
C...          ........                                  ........
C...          left bulb                                right bulb
C...
C...  We make the following assumptions:
C...
C...  (1)  The pressure and temperature in the system are uniform.  Also,
C...       the ideal gas law can be used to calculate the total molar
C...       density.
C...
C...  (2)  The diffusion in the tube is essentially at steady state,
C...       i.e., the concentration profiles in the tube adjust in-
C...       stantaneously to changes in the concentrations in the
C...       bulbs.
C...
C...  (3)  The concentration profiles in the tube are linear.
C...
C...  (4)  The fluxes of the three components in the tubes are evaluated
C...       at the midpoint.
C...
C...  (5)  The concentrations in the bulbs are uniform (i.e., dependent
C...       on time only, and not position).
C...
C...  The following three component system is considered:
C...
C...     i = 1; hydrogen
C...
C...     i = 2; nitrogen
C...
C...     i = 3; carbon dioxide
C...
C...  The objective then is to compute the concentrations in the bulbs
C...  as a function of time.
C...
C...  The following equations are used to model the system (capital
C...  letters are used for consistency with the Fortran programming)
C...
C...     Component balances on the left bulb:
C...
C...        dXIO/dt = -NI*A/(CT*VO), I = 1,2                            (1)
C...
C...  where
C...
C...     XIO   mole fraction of component i (= I) in the left bulb (gm
C...           mols I/gm mol gas)
C...
C...     t     time (s)
C...
C...     NI    flux of component I in the tube (gm mols of I/s-cm**2)
C...
C...     A     tube cross sectional area (cm**2)
C...
C...     CT    total molar density (gm mols gas/cm**3)
C...
C...     VO    volume of the left bulb (cm**3)
C...
C...  Similarly a component balance on the right bulb gives (with "O" in
C...  eq. (1) replaced by "L")
C...
C...        dXIL/dt = NL*A/(CT*VL), I = 1,2                             (2)
C...
C...  The third component concentration and flux can be computed from
C...
C...     X1O + X2O + X3O = 1                                            (3)
C...
C...     X1L + X2L + X3L = 1                                            (4)
C...
C...     N1 + N2 + N3 = 0                                               (5)
C...
C...  The fluxes are related to the gradients by the Stefan-Maxwell
C...  equation [Bird, et al (1960), pp 569-572]
C...
C...                NC XI*NJ - XJ*NI
C...     DEL(XI) = SUM -------------                                    (6)
C...               J=1    CT*DIJ
C...              I NE J
C...  where
C...
C...     DEL   "del" differential operator
C...
C...     NC    number of components (= 3)
C...
C...     DIJ   binary diffusion coefficient for components I and J
C...
C...  For I = 1, eq. (6) is
C...
C...                3  XI*NJ - XJ*NI
C...     DEL(X1) = SUM -------------
C...               J=2     CT*DIJ
C...              J NE 1
C...
C...     = (1/(CT*D12))*(X1*N2 - X2*N1) + (1/(CT*D13))*(X1*N3 - X3*N1)  (7)
C...
C...  Eq. (7) can be written in standard algebraic form as
C...
C...     A11*N1 + A12*N2 = B1                                           (8)
C...
C...  where
C...
C...     A11 = (1/(CT*D12))*(-X2) + (1/(CT*D13))*(-X1-X3)
C...
C...     A12 = (1/(CT*D12))*( X1) + (1/(CT*D13))*(-X1)
C...
C...      B1 = DEL(X1)
C...
C...  For I = 2, eq. (6) is
C...
C...                3  XI*NJ - XJ*NI
C...     DEL(X2) = SUM -------------
C...               J=1     CT*DIJ
C...              J NE 2
C...
C...     = (1/(CT*D21))*(X2*N1 - X1*N2) + (1/(CT*D23))*(X2*N3 - X3*N2)  (9)
C...
C...  Eq. (9) can be written in standard algebraic form as
C...
C...     A21*N1 + A22*N2 = B2                                          (10)
C...
C...  where
C...
C...     A21 = (1/(CT*D21))*( X2) + (1/(CT*D23))*(-X2)
C...
C...     A22 = (1/(CT*D21))*(-X1) + (1/(CT*D23))*(-X2-X3)
C...
C...      B2 = DEL(X2)
C...
C...  For the H2-N2-CO2 SYSTEM, the following diffusivities (cm**2/s)
C...  are reported by Taylor, et al (1993), pp 105-110, 131-133
C...
C...     D12 = 8.33D-01
C...
C...     D13 =  6.8D-01
C...
C...     D23 = 1.68D-01
C...
C...  Eqs. (1) to (11) are programmed in the following coding to compute
C...  the three component concentrations and fluxes.
C...
C...  Bird, R. B., (1960), Transport Phenomena, Wiley, New York
C...
C...  Taylor, R., and R. Krishna, (1993), Multicomponent Mass Transfer,
C...  Wiley, New York
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NOSTP,  NORUN
     +      /Y/    X1O,    X2O,    X1L,    X2L
     +      /F/   X1OT,   X2OT,   X1LT,   X2LT
     +      /C/      A,     ZL,     VO,     VL,
     +              TC,      P,      R,     CT,
     +     D12,    D13,    D23,     X1,     X2,
     +             X3O,    X3L,  DELX1,  DELX2,
     +              N1,     N2,     N3
     +      /I/     IP
      DOUBLE PRECISION      N1,     N2,     N3
C...
C...  MODEL PARAMETERS
C...
C...     TUBE LENGTH (CM)
         ZL=10.0D0
C...
C...     TUBE CROSS SECTIONAL AREA (CM**2)
         A=0.04D0
C...
C...     TERMINAL VOLUMES (CM**3)
         VO=80.0D0
         VL=80.0D0
C...
C...     TEMPERATURE (C)
         TC=35.2D0
C...
C...     PRESSURE (PA)
         P=101300.0D0
C...
C...     GAS CONSTANT (PA-CM**3/K-GM MOL)
         R=101300.0D0*22400.0D0/273.16D0
C...
C...     MOLAR DENSITY (GM MOL/CM**3)
         CT=P/(R*(TC+273.16D0))
C...
C...     BINARY DIFFUSIVITIES (CM**2/HR)
         D12=8.33D-01*3600.0D0
         D13= 6.8D-01*3600.0D0
         D23=1.68D-01*3600.0D0
C...
C...  INITIAL CONDITIONS
      X1O=0.0D0
      X2O=0.5D0
      X1L=0.5D0
      X2L=0.5D0
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/      T,  NOSTP,  NORUN
     +      /Y/    X1O,    X2O,    X1L,    X2L
     +      /F/   X1OT,   X2OT,   X1LT,   X2LT
     +      /C/      A,     ZL,     VO,     VL,
     +              TC,      P,      R,     CT,
     +     D12,    D13,    D23,     X1,     X2,
     +             X3O,    X3L,  DELX1,  DELX2,
     +              N1,     N2,     N3
     +      /I/     IP
      DOUBLE PRECISION      N1,     N2,     N3
C...
C...  GRADIENTS
      DELX1=(X1L-X1O)/ZL
      DELX2=(X2L-X2O)/ZL
C...
C...  MIDPOINT CONCENTRATIONS
      X1=(X1L+X1O)/2.0D0
      X2=(X2L+X2O)/2.0D0
C...
C...  DEFINE LINEAR ALGEBRAIC EQUATIONS FOR N1 AND N2
      X3=1.0D0-X1-X2
      A11=1.0D0/(CT*D12)*(-X2)+(1.0D0/(CT*D13))*(-X1-X3)
      A12=1.0D0/(CT*D12)*( X1)+(1.0D0/(CT*D13))*(-X1)
       B1=DELX1
      A21=1.0D0/(CT*D12)*( X2)+(1.0D0/(CT*D23))*(-X2)
      A22=1.0D0/(CT*D12)*(-X1)+(1.0D0/(CT*D23))*(-X2-X3)
       B2=DELX2
C...
C...  SOLUTION OF LINEAR ALGEBRAIC EQUATIONS BY KRAMER'S RULE
C...
C...     DISCRIMINANT
         DIS=A11*A22-A12*A21
C...
C...     CHECK FOR NEAR-SINGULAR SYSTEM
         IF(DABS(DIS).LT.1.0D-20)THEN
            WRITE(*,1)DIS,T
1           FORMAT(' Discriminant = ',D12.3,' at t = ',F6.2)
            STOP
         END IF
C...
C...     FLUXES
         N1=(A22*B1-A12*B2)/DIS
         N2=(A11*B2-A21*B1)/DIS
C...
C...  ODES
      X1OT=-N1*A/(CT*VO)
      X2OT=-N2*A/(CT*VO)
      X1LT= N1*A/(CT*VL)
      X2LT= N2*A/(CT*VL)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NT=101)
      COMMON/T/      T,  NOSTP,  NORUN
     +      /Y/    X1O,    X2O,    X1L,    X2L
     +      /F/   X1OT,   X2OT,   X1LT,   X2LT
     +      /C/      A,     ZL,     VO,     VL,
     +              TC,      P,      R,     CT,
     +     D12,    D13,    D23,     X1,     X2,
     +             X3O,    X3L,  DELX1,  DELX2,
     +              N1,     N2,     N3
     +      /I/     IP
      DOUBLE PRECISION      N1,     N2,     N3
C...
C...  MONITOR THE OUTPUT
      IP=IP+1
      WRITE(*,*)IP,T
C...
C...  HEADING FOR THE SOLUTION
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(9X,'t',7X,'X1O',7X,'X2O',7X,'X3O',/,
     +             17X,'X1L',7X,'X2L',7X,'X3L',/,
     +             17X,' N1',7X,' N2',7X,' N3',/,
     + 3X,'dX1O/dt',3X,'dX2O/dt',3X,'dX1L/dt',3X,'dX2L/dt')
C...
C...  WRITE THE SOLUTION EVERY 10TH CALL TO PRINT
      X3O=1.0D0-X1O-X2O
      X3L=1.0D0-X1L-X2L
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
         WRITE(NO,2)T,X1O,X2O,X3O,X1L,X2L,X3L,
     +              N1,N2,N3,X1OT,X2OT,X1LT,X2LT
2        FORMAT(F10.2,3F10.5,/,
     +            10X,3F10.5,/,
     +            10X,3D10.2,/,
     +                4D10.2,//)
      END IF
C...
C...  PLOT THE SOLUTION WITH MATLAB
      CALL MPLOT
      RETURN
      END
      SUBROUTINE MPLOT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NT=101)
      COMMON/T/      T,  NOSTP,  NORUN
     +      /Y/    X1O,    X2O,    X1L,    X2L
     +      /F/   X1OT,   X2OT,   X1LT,   X2LT
     +      /C/      A,     ZL,     VO,     VL,
     +              TC,      P,      R,     CT,
     +     D12,    D13,    D23,     X1,     X2,
     +             X3O,    X3L,  DELX1,  DELX2,
     +              N1,     N2,     N3
     +      /I/     IP
      DOUBLE PRECISION      N1,     N2,     N3
C...
C...  OPEN FILES FOR MATLAB PLOTTING
      IF(IP.EQ.1)THEN
         OPEN(1,FILE='ma3c.out')
         OPEN(2,FILE='ma3f.out')
      END IF
C...
C...  WRITE THE SOLUTION FOR PLOTTING
C...
C...     CONCENTRATIONS
         X3O=1.0D0-X1O-X2O
         X3L=1.0D0-X1L-X2L
         WRITE(1,1)T, X1O, X2O, X3O, X1L, X2L, X3L
1        FORMAT(7F9.5)
C...
C...     FLUXES
         N3=-N1-N2
         WRITE(2,2)T, N1, N2, N3
2        FORMAT(F8.4,3F13.7)
      RETURN
      END
