*DECK SEP07
      SUBROUTINE INITAL
C...
C...  MULTICOMPONENT CHROMATOGRAPHIC SEPARATION
C...
C...  THE FOLLOWING PARTIAL DIFFERENTIAL EQUATIONS (PDES) ARE DERIVED
C...  FOR AN ISOTHERMAL CHROMATOGRAPHIC COLUMN WHICH IS INTENDED FOR
C...  THE SEPARATION OF TWO COMPONENTS INJECTED AT THE ENTRANCE OF THE
C...  COLUMN
C...
C...  COMPONENT 1 FLUID BALANCE
C...
C...     A*DZ*EPS*RHOF*C1  = -V*A*EPS*RHOF*(C1     - C1 )
C...                     T                    Z+DZ     Z
C...                                                                (1)
C...                                                         E
C...                         -A*DZ*(1-EPS)*(RHOS/MW1)*KP1*(Q1  - Q1)
C...
C...  THIS EQUATION IS A COMPONENT 1 FLUID BALANCE ON A SECTION OF THE
C...  COLUMN OF LENGTH DZ.  THE TERMS IN THE EQUATION ARE
C...
C...     A*DZ*EPS*RHOF*C1
C...                     T
C...
C...        ACCUMULATION OF COMPONENT 1 IN THE FLUID IN THE SECTION OF
C...        LENGTH DZ
C...
C...        UNITS     (CM**2 COLUMN)*(M COLUMN)*(CM**3 FLUID/CM**3 COLUM
C...                  (MOLS FLUID/CM**3 FLUID)*(MOLS COMP 1/MOL FLUID)*
C...                  (1/SEC) = MOLS COMP 1/SEC
C...
C...     -V*A*EPS*RHOF*(C1     - C1 )
C...                      Z+DZ     Z
C...
C...        NET FLOW OF COMPONENT 1 INTO THE SECTION OF LENGTH DZ
C...
C...        UNITS     (M COLUMN/SEC)*(CM**2 COLUMN)*(CM**3 FLUID/CM**3
C...                  COLUMN)*(MOLS FLUID/CM**3 FLUID)*(MOLS COMP 1/
C...                  MOLS FLUID) = MOLS COMP 1/SEC
C...
C...                                     E
C...     -A*DZ*(1-EPS)*(RHOS/MW1)*KP1*(Q1  - Q1)
C...
C...        NET TRANSFER OF COMPONENT 1 FROM THE SOLID TO THE FLUID
C...
C...        UNITS     (CM**2 COLUMN)*(M COLUMN)*(CM**3 SOLID/CM**3 COLUM
C...                  *(GM SOLID/CM**3 SOLID)*(MOLS COMP 1/GM COMP 1)
C...                  *(GM COMP 1/GM SOLID) = MOLS COMP 1/SEC
C...
C...  WHERE
C...
C...     C1     MOLAR CONCENTRATION OF COMPONENT 1 (MOLS COMP 1/MOL
C...            FLUID)
C...
C...     Z      AXIAL POSITION ALONG THE COLUMN (M COLUMN)
C...
C...     T      TIME (SEC)
C...
C...     A      COLUMN CROSS SECTIONAL AREA (CM**2 COLUMN)
C...
C...     EPS    COLUMN VOID FRACTION (CM**3 FLUID/CM**3 COLUMN)
C...
C...     RHOF   FLUID MOLAR DENSITY (MOLS FLUID/CM**3 FLUID)
C...
C...     RHOS   SOLID MASS DENSITY (GM SOLID/CM**3 SOLID)
C...
C...     V      FLUID VELOCITY (M COLUMN/SEC)
C...
C...     KP1    FLUID TO SOLID MASS TRANSFER COEFFICIENT FOR COMPONENT 1
C...            (1/SEC)
C...
C...     MW1    MOLECULAR WEIGHT OF COMPONENT 1 (GM COMP 1/MOL COMP 1)
C...
C...       E    EQUILIBRIUM SOLID CONCENTRATION OF COMPONENT 1 (GM COMP
C...     Q1     1/GM SOLID)
C...
C...     Q1     SOLID CONCENTRATION OF COMPONENT 1 (GM COMP 1/GM SOLID)
C...
C...  DIVISION BY A*DZ*EPS*RHOF, FOLLOWED BY DZ ---+ 0 GIVES
C...
C...                                                       E
C...     C1  = -V*C1  - (1-EPS)*RHOS/(EPS*RHOF*MW1)*KP1*(Q1  - Q1)  (2)
C...       T        Z
C...
C...  COMPONENT 1 SOLID BALANCE
C...
C...                                                     E
C...    A*DZ*(1-EPS)*RHOS*Q1  = A*DZ*(1-EPS)*RHOS*KP1*(Q1  - Q1)    (3)
C...                        T
C...
C...  THIS EQUATION IS A COMPONENT 1 SOLID BALANCE ON A SECTION OF THE
C...  COLUMN OF LENGTH DZ.  THE TERMS IN THE EQUATION ARE
C...
C...     A*DZ*(1-EPS)*RHOS*Q1
C...                         T
C...
C...        ACCUMULATION OF COMPONENT 1 ON THE SOLID IN THE SECTION OF
C...        LENGTH DZ
C...
C...        UNITS    (CM**2 COLUMN)*(M COLUMN)*(CM**3 SOLID/CM**3 COLUMN
C...                 (GM SOLID/CM**3 SOLID)*(GM COMP 1/GM SOLID)*(1/SEC)
C...                 = GM COMP 1/SEC
C...
C...                              E
C...     A*DZ*(1-EPS)*RHOS*KP1*(Q1  - Q1)
C...
C...        NET FLOW OF COMPONENT 1 FROM THE SOLID TO THE FLUID
C...
C...        UNITS    (CM**2 COLUMN)*(M COLUMN)*(CM**3 SOLID/CM**3 COLUMN
C...                 (GM SOLID/CM**3 SOLID)*(1/SEC)*(GM COMP 1/GM SOLID)
C...                 = GM COMP 1/SEC
C...
C...  WHERE ALL OF THE PARAMETERS AND VARIABLES WERE DEFINED PREVIOUSLY.
C...
C...  DIVISION BY A*DZ*(1-EPS)*RHOS GIVES
C...
C...                  E
C...     Q1  = KP1*(Q1  - Q1)                                       (4)
C...       T
C...
C...  EQUATION (2) REQUIRES AN INITIAL CONDITION AND A BOUNDARY
C...  CONDITION
C...
C...     C1(0,T) = F1(T), C1(Z,0) = C10                          (5)(6)
C...
C...  WHILE EQUATION (4) REQUIRES AN INITIAL CONDITION
C...
C...     Q1(Z,0) = Q10    (7)
C...
C...  ALSO, THE EQUILIBRIUM SOLID CONCENTRATION IS GIVEN BY A LINEAR
C...  ISOTHERM
C...
C...       E
C...     Q1  = QA1*K1*C1                                            (8)
C...
C...  A SIMILAR SET OF EQUATIONS CAN BE WRITTEN FOR COMPONENT 2.  THEN
C...  THE MODEL FOR THE CHROMATOGRAPHIC SYSTEM CONSISTS OF TWO PARTIAL
C...  DIFFERENTIAL EQUATIONS (PDES) FOR THE FLUID, TWO PDES FOR THE
C...  SOLID, TWO BOUNDARY CONDITIONS, FOUR INITIAL CONDITIONS AND TWO
C...  EQUILIBRIUM RELATIONSHIPS.
C...
C...  THE MODEL PARAMETERS ARE
C...
C...     EPS = 0.360, ZL = 100, V = 1.0, RHOF = 0.01, RHOS = 0.8,
C...
C...     KP1 = KP2 = 1, K1 = 1, K2 = 0.3, QA1 = 2.90, QA2 = 3.67
C...
C...     C10 = C20 = 1.0E-15, Q10 = Q20 = 1.0E-15, MW1 = MW2 = 142
C...
C...     N = 21
C...
C...     F1(T) = F2(T) = 0.03, 0 LE T LE 30
C...
C...     F1(T) = F2(T) = 0, T GT 30
C...
C...  THE EXITING CONCENTRATIONS, C1(ZL,T) AND C2(ZL,T), ARE OF PARTI-
C...  CULAR INTEREST IN THIS SIMULATION.
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    C1(21),    C2(21),    Q1(21),    Q2(21)
     2      /F/ DC1DT(21), DC2DT(21), DQ1DT(21), DQ2DT(21)
     3      /S/ DC1DZ(21), DC2DZ(21),   QE1(21),   QE2(21)
     4      /C/         N,        ZL,       QA1,       QA2,
     5                 P1,        P2,        P3,        P4,
     6                 K1,        K2,         V,        IP
      REAL            KP1,       KP2,        K1,        K2,
     1                MW1,       MW2
C...
C...  SET THE MODEL PARAMETERS
      EPS=0.360
      ZL=100.
      V=1.
      RHOF=0.01
      RHOS=0.8
      KP1=1.0
      KP2=1.0
      K1=1.0
      K2=0.3
      QA1=2.90
      QA2=3.67
      MW1=142.
      MW2=142.
      C10=1.0E-15
      C20=1.0E-15
      Q10=1.0E-15
      Q20=1.0E-15
      N=21
      IP=0
C...
C...  COMPUTE SOME CONSTANTS TO BE USED IN DERV
      P1=(1.-EPS)*RHOS/(EPS*RHOF*MW1)*KP1
      P2=(1.-EPS)*RHOS/(EPS*RHOF*MW2)*KP2
      P3=KP1
      P4=KP2
      QA1=QA1*K1
      QA2=QA2*K2
C...
C...  SET THE INITIAL FLUID CONCENTRATIONS
      DO 1 I=1,N
      C1(I)=C10
      C2(I)=C20
      Q1(I)=Q10
      Q2(I)=Q20
1     CONTINUE
C...
C...  INITIALIZE THE DERIVATIVE CALCULATIONS
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    C1(21),    C2(21),    Q1(21),    Q2(21)
     2      /F/ DC1DT(21), DC2DT(21), DQ1DT(21), DQ2DT(21)
     3      /S/ DC1DZ(21), DC2DZ(21),   QE1(21),   QE2(21)
     4      /C/         N,        ZL,       QA1,       QA2,
     5                 P1,        P2,        P3,        P4,
     6                 K1,        K2,         V,        IP
      REAL            KP1,       KP2,        K1,        K2
C...
C...  MODEL BOUNDARY CONDITIONS
      C1(1)=0.
      C2(1)=0.
      IF(T.LE.30.)C1(1)=0.03
      IF(T.LE.30.)C2(1)=0.03
C...
C...  SET CONSTRAINTS ON THE CONCENTRATIONS
      DO 1 I=1,N
      IF(C1(I).LT.0.)C1(I)=0.
      IF(C2(I).LT.0.)C2(I)=0.
      IF(Q1(I).LT.0.)Q1(I)=0.
      IF(Q2(I).LT.0.)Q2(I)=0.
1     CONTINUE
C...
C...  COMPUTE THE EQUILIBRIUM SOLID CONCENTRATIONS
      CALL EQUIL(N,C1,QE1,QA1,K1)
      CALL EQUIL(N,C2,QE2,QA2,K2)
C...
C...  COMPUTE THE SPATIAL DERIVATIVES
      CALL DSS020(0.,ZL,N,C1,DC1DZ,1.)
      CALL DSS020(0.,ZL,N,C2,DC2DZ,1.)
C...
C...  MODEL PDES
      DO 2 I=1,N
      DC1DT(I)=-V*DC1DZ(I)-P1*(QE1(I)-Q1(I))
      DC2DT(I)=-V*DC2DZ(I)-P2*(QE2(I)-Q2(I))
      DQ1DT(I)=P3*(QE1(I)-Q1(I))
      DQ2DT(I)=P4*(QE2(I)-Q2(I))
2     CONTINUE
      DC1DT(1)=0.
      DC2DT(1)=0.
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    C1(21),    C2(21),    Q1(21),    Q2(21)
     2      /F/ DC1DT(21), DC2DT(21), DQ1DT(21), DQ2DT(21)
     3      /S/ DC1DZ(21), DC2DZ(21),   QE1(21),   QE2(21)
     4      /C/         N,        ZL,       QA1,       QA2,
     5                 P1,        P2,        P3,        P4,
     6                 K1,        K2,         V,        IP
      REAL            KP1,       KP2,        K1,        K2
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION   TP(101), CP(4,101)
C...
C...  PRINT THE SOLUTION AT THE BEGINNING AND END OF EACH RUN
      IF((IP.NE.0).AND.(IP.NE.100))GO TO 2
      N2=1
      WRITE(NO,1)T,(C1(I),I=1,N,N2),(C2(I),I=1,N,N2),
     1             (Q1(I),I=1,N,N2),(Q2(I),I=1,N,N2)
1     FORMAT(11H          T,F7.1,//,
     1       11H    C1(Z,T),E11.3,/,5E11.3,/,5E11.3,/,
     2                              5E11.3,/,5E11.3,//,
     3       11H    C2(Z,T),E11.3,/,5E11.3,/,5E11.3,/,
     4                              5E11.3,/,5E11.3,//,
     5       11H    Q1(Z,T),E11.3,/,5E11.3,/,5E11.3,/,
     6                              5E11.3,/,5E11.3,//,
     7       11H    Q2(Z,T),E11.3,/,5E11.3,/,5E11.3,/,
     8                              5E11.3,/,5E11.3,//)
C...
C...  STORE THE NUMERICAL SOLUTION FOR PLOTTING
2     IP=IP+1
      CP(1,IP)=C1(1)
      CP(2,IP)=C2(1)
      CP(3,IP)=C1(N)
      CP(4,IP)=C2(N)
      TP(IP)=T
C...
C...  PLOT THE NUMERICAL SOLUTION
      IF(IP.LT.101)RETURN
      CALL TPLOTS(4,IP,TP,CP)
      WRITE(NO,3)
3     FORMAT(1H ,//,13H 1 - C1(0.,T),3X,13H 2 - C2(0.,T),/,
     1              13H 3 - C1(ZL,T),3X,13H 4 - C2(ZL,T))
5     CONTINUE
C...
C...  ******************************************************************
C...
C...  THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE-DEPENDENT
C...  AND WILL HAVE TO BE CHANGED FOR THE USER*S LOCAL COMPUTER OR DE-
C...  ACTIVATED (CONVERTED TO COMMENTS)
C...
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,0.,100.,7.,-0.005,0.005)
C...  CALL QIKPLT(TP,CP,101,12H*TIME (SEC)*,
C... 1            29H*C1(Z,T), C2(Z,T), Z = 0, ZL*,2H**,4)
C...
C...  ******************************************************************
C...
      RETURN
      END
      SUBROUTINE EQUIL(N,C,QE,QA,K)
C...
C...  SUBROUTINE COMPUTES THE EQUILIBRIUM SOLID CONCENTRATION
C...
C...  ARGUMENT LIST
C...
C...     N        NUMBER OF COMPONENTS (INPUT)
C...
C...     C        ONE-DIMENSIONAL ARRAY OF FLUID CONCENTRATIONS
C...              (INPUT)
C...
C...     QE       ONE-DIMENSIONAL ARRAY OF EQUILIBRIUM SOLID
C...              CONCENTRATIONS (OUTPUT)
C...
C...     QA       SOLID CAPACITANCE (INPUT)
C...
C...     K        CONSTANT IN THE ADSORPTION ISOTHERM (INPUT)
C...

      DIMENSION C(N),QE(N)
      REAL K
C...
C...  COMPUTE THE EQUILIBRIUM SOLID CONCENTRATIONS FOR THE N COMPONENTS
C...  (NOTE - A SIMPLE LINEAR EQUILIBRIUM ISOTHERM IS USED IN DO LOOP
C...  1.  HOWEVER, ANY ADSORPTION ISOTHERM CAN BE PROGRAMMED, INCLUDING
C...  A MULTICOMPONENT, NONLINEAR ISOTHERM.  PERHAPS ADDITIONAL ADSORP-
C...  TION CONSTANTS WILL HAVE TO BE PASSED TO THE SUBROUTINE FROM THE
C...  CALLING PROGRAM, BUT THE COMPLEXITY OF THE ISOTHERM IS, IN PRIN-
C...  CIPLE, ESSENTIALLY UNLIMITED)
      DO 1 I=1,N
      QE(I)=QA*C(I)
1     CONTINUE
      RETURN
      END

