*DECK SEP35
      SUBROUTINE INITAL
C...
C...  DYNAMIC SIMULATION OF A PACKED COLUMN ADSORBER
C...
C...  Develop a dynamic simulation for a packed colum adsorber using
C...  the following variables:
C...
C...  y(z,t)      gas phase mole fraction of component to be adsorbed
C...              (gm mol adsorbate/mol gas)
C...
C...  q(z,t)      solid phase concentration of component to be adsorbed
C...              (gm mol adsorbate/gm solid)
C...
C...  g           molar flow rate of gas through the adsorber (gm mol
C...              gas/sec)
C...
C...  eps         fraction of solid adsorbent in the column (cm**3
C...              solid/cm**3 column)
C...
C...  rhog        molar density of the gas (gm mol gas/cm**3)
C...
C...  rhos        density of the solid adsorbent (gm solid/cm**3 solid)
C...
C...  A           cross sectional area of the column (cm**2 column)
C...
C...  z           axial position along the adsorber (cm column)
C...
C...  t           time (sec)
C...
C...  zl          length of the adsorber (cm column)
C...
C...  The rate of mass transfer from the gas to the solid phase is given
C...  by
C...
C...     r = c1*y*(qe - q) - c2*q  (gm mol adsorbate/gm solid-sec)
C...
C...  where c1 and c2 are given rate constants for the solid adsorption
C...  and qe is the solid concentration in equilibrium with the gas of
C...  composition y given by the nonlinear adsorption isotherm:
C...
C...     qe = c3*y/(1 + c4*y) (gm mol adsorbate/gm solid)
C...
C...  c3 and c4 are given adsorption isotherm constants.
C...
C...  The model should consist of the partial differential equations
C...  (PDEs) for y(z,t) and q(z,t) with the appropriate initial and
C...  boundary conditions.  Initially, the column starts out with no
C...  gas to be adsorbed in either the gas or solid phases.
C...
C...  A gas-phase balance on a section of the column of length dz gives
C...
C...     (1 - eps)*A*dz*rhog*y  = g*y  - g*y     - A*dz*eps*rhos*r
C...                          t      z      z+dz
C...
C...  Division by the coefficient of the time derivative, followed by
C...  dz ---+ 0 gives
C...
C...     y  = -v*y  - (eps*rhos/((1 - eps)*rhog))*r                (1)
C...      t       z
C...
C...     v = g/((1 - eps)*A*rhog) (A = column cross sectional area)
C...
C...  Similarly, a balance on the solid phase gives
C...
C...     A*dz*eps*rhos*q  = A*dz*eps*rhos*r
C...                    t
C...
C...  or, after dividing by the coefficient of the time derivative
C...
C...     q  = r                                                    (2)
C...      t
C...
C...  The initial and boundary conditions for equations (1) and (2) are
C...
C...     y(z,0) = y0, y(0,t) = yin(t), q(z,0) = q0           (3)(4)(5)
C...
C...  Equations (1) to (5) are programmed in the following coding:
C...
C...  DSS/2 common area
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     Y(21),     Q(21)
     2      /F/    YT(21),    QT(21)
     3      /S/    YZ(21)
     4      /C/         V,        C1,        C2,        C3,        C4,
     5                 C5,       YIN,    QE(21),     R(21),        ZL,
     6                  N,        IP
C...
C...  SET THE MODEL PARAMETERS
      A=100.
      ZL=100.
      RHOG=0.1
      RHOS=1.0
      EPS=0.6
      C1=100.
      C2=10.
      C3=100.
      C4=1000.
      G=30.
      V=G/((1.0-EPS)*A*RHOG)
      C5=EPS*RHOS/((1.0-EPS)*RHOG)
      YIN=0.01
C...
C...  INITIAL CONDITIONS (3) AND (5)
      N=21
      Y0=0.
      Q0=0.
      DO 1 I=1,N
      Y(I)=Y0
      Q(I)=Q0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     Y(21),     Q(21)
     2      /F/    YT(21),    QT(21)
     3      /S/    YZ(21)
     4      /C/         V,        C1,        C2,        C3,        C4,
     5                 C5,       YIN,    QE(21),     R(21),        ZL,
     6                  N,        IP
C...
C...  LIMIT THE CONCENTRATIONS TO NON-NEGATIVE VALUES
      CALL LIMITS
C...
C...  BOUNDARY CONDITION (4)
      Y(1)=YIN
C...
C...  DERIVATIVE Y
C...              Z
      CALL DSS020(0.,ZL,N,Y,YZ,1.0)
C...
C...  EQUATIONS (1) AND (2)
      DO 1 I=1,N
C...
C...     EQUILIBRIUM SOLID CONCENTRATION
         QE(I)=C3*Y(I)/(1.0+C4*Y(I))
C...
C...     RATE OF ADSORPTION
         R(I)=C1*Y(I)*(QE(I)-Q(I))-C2*Q(I)
C...
C...     YT AT Z = 0
            IF(I.EQ.1)THEN
               YT(1)=0.
C...
C...     YT AT Z NE 0
            ELSE IF (I.NE.1)THEN
               YT(I)=-V*YZ(I)-C5*R(I)
            END IF
C...
C...     QT
            QT(I)=R(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE LIMITS
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     Y(21),     Q(21)
     2      /F/    YT(21),    QT(21)
     3      /S/    YZ(21)
     4      /C/         V,        C1,        C2,        C3,        C4,
     5                 C5,       YIN,    QE(21),     R(21),        ZL,
     6                  N,        IP
C...
C...  LIMIT THE CONCENTRATIONS TO NON-NEGATIVE VALUES
      DO 1 I=1,N
      IF( Y(I).LT.0.) Y(I)=0.
      IF( Q(I).LT.0.) Q(I)=0.
      IF(QE(I).LT.0.)QE(I)=0.
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     Y(21),     Q(21)
     2      /F/    YT(21),    QT(21)
     3      /S/    YZ(21)
     4      /C/         V,        C1,        C2,        C3,        C4,
     5                 C5,       YIN,    QE(21),     R(21),        ZL,
     6                  N,        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      REAL TP(101), YP(2,101)
C...
C...  PRINT THE SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT
      IP=IP+1
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
         WRITE(NO,1)T,( Y(I),I=1,N,1),( Q(I),I=1,N,1),(QE(I),I=1,N,1),
     1                ( R(I),I=1,N,1),(YT(I),I=1,N,1),(QT(I),I=1,N,1)
1     FORMAT(1H ,////,' T = ',F5.1,//,
     1 '   Y',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//,
     1 '   Q',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//,
     1 '  QE',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//,
     1 '   R',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//,
     1 '  YT',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//,
     1 '  QT',5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,5E10.3,/,4X,E10.3,//)
C...
C...  STORE THE SOLUTION FOR PLOTTING
      TP(IP)=T
      YP(1,IP)=Y(1)
      YP(2,IP)=Y(N)
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.101)RETURN
      CALL TPLOTS(2,IP,TP,YP)
      WRITE(NO,11)
11    FORMAT(1H ,//,' 1 - Y(0,T)   2 - Y(ZL,T) VS T')
      END IF
      RETURN
      END

