      SUBROUTINE INITAL
C...
C...  DIFFUSION AND NONLINEAR ADSORPTION IN A PORE
C...
C...  CONSIDER DIFFUSION INTO AND ADSORPTION ONTO AN IDEALIZED PORE
C...
C             .............................................
C                                                         .
C     C(0,T) = CS                  C(Z,T)                 .C (L,T) = 0
C                                  G(Z,T)                 . Z
C             .............................................
C           Z = 0                                       Z = L
C             +..................... L ...................+
C
C...  A MATERIAL BALANCE ON THE FLUID PHASE IN A DIFFERENTIAL LENGTH OF
C...  THE PORE GIVES
C...
C...          C  = D*C  - (2/(R0*SM))*(KA*(1 - G)*C - KD*G)          (1)
C...           T      ZZ
C...
C...  SIMILARLY, A MATERIAL BALANCE ON A DIFFERENTIAL LENGTH OF THE PORE
C...  WALL GIVES
C...
C...          G  = KA*(1 - G)*C - KD*G                               (2)
C...           T
C...
C...  WHERE
C...
C...     C(Z,T)    CONCENTRATION OF ADSORBATE IN THE FLUID AT POSITION Z
C...               ALONG THE PORE AND TIME T
C...
C...     G(Z,T)    FRACTION OF THE PORE WALL COVERED AT POSITION Z ALONG
C...               THE PORE AND TIME T
C...
C...     R0        RADIUS OF THE PORE
C...
C...     L         LENGTH OF THE PORE
C...
C...     SM        SURFACE AREA OF THE PORE OCCUPIED BY ONE MOLE OF
C...               ADSORBATE
C...
C...     D         DIFFUSIVITY OF THE ADSORBATE IN THE FLUID WITHIN THE
C...               PORE
C...
C...     KA,KD     KINETIC CONSTANTS FOR ADSORPTION AND DESORPTION,
C...               RESPECTIVELY
C...
C...     CS        CONCENTRATION OF ADSORBATE AT THE MOUTH (Z = 0) OF
C...               THE PORE
C...
C...  THE BOUNDARY CONDITIONS FOR EQUATION (1) REFLECT THE CONCENTRATION
C...  CS AT THE MOUTH OF THE PORE AND NO DIFFUSION THROUGH THE END OF
C...  THE PORE
C...
C...          C(0,T) = CS, C (L,T) = 0                            (3)(4)
C...                        Z
C...
C...  INITIALLY, THE FLUID WITHIN THE PORE AND THE PORE WALL ARE FREE OF
C...  ADSORBATE
C...
C...          C(Z,0) = 0, G(Z,0) = 0                              (5)(6)
C...
C...  EQUATIONS (1) TO (6) CONSTITUTE THE PROBLEM TO BE SOLVED, I.E.,
C...  C(Z,T) AND G(Z,T) ARE TO BE COMPUTED NUMERICALLY.
C...
C...  TO FACILITATE THE SOLUTION, DIMENSIONLESS VARIABLES ARE DEFINED AS
C...
C...          F = C/CS, Z = Z/L, T = T*D/L**2
C...
C...  NOTE THAT T NOW REPRESENTS DIMENSIONLESS TIME.  EQUATIONS (1) TO
C...  (6) CAN BE STATED IN TERMS OF F, Z AND T AS
C...
C...          F  = F   + HD*G - HA*(1 - G)*F                         (7)
C...           T    ZZ
C...
C...          G  = H*(HA*(1 - G)*F - HD*G)                           (8)
C...           T
C...
C...          F(0,T) = 1, F (1,T) = 0                            (9)(10)
C...                       Z
C...
C...          F(Z,0) = 0, G(Z,0) = 0                            (11)(12)
C...
C...  WHERE
C...
C...     HA        2*(L**2)*KA/(D*R0*SM)
C...
C...     HD        2*(L**2)*KD/(CS*D*R0*SM)
C...
C...     H         CS*R0*SM/2
C...
C...  THE METHOD OF LINES IS IMPLEMENTED VIA THREE APPROACHES (NCASE =
C...  1, 2, 3).  THE DERIVATIVE F   IN EQUATION (7) IS COMPUTED VIA
C...                             ZZ
C...
C...     (1)  A THREE-POINT CENTRAL DIFFERENCE APPROXIMATION EXPLICITLY
C...          PROGRAMMED IN SUBROUTINE DERV1
C...
C...     (2)  THE THREE-POINT CENTRAL DIFFERENCE APPROXIMATION IN
C...          SUBROUTINE DSS002 CALLED IN DERV2
C...
C...     (3)  THE FIVE-POINT CENTRAL DIFFERENCE APPROXIMATION IN
C...          SUBROUTINE DSS004 CALLED IN DERV3
C...
C...  AN INTERESTING INFINITE-TIME (STEADY STATE) SOLUTION TO EQUATIONS
C...  (7) TO (12) CAN ALSO BE DERIVED AS A CHECK ON THE NUMERICAL SOLU-
C...  TION.  THUS, FOR F  = 0, G  = 0, EQUATION (8) CAN BE SOLVED
C...                    T       T
C...  EXPLICITLY FOR G (NOTE THAT THIS IS POSSIBLE BECAUSE EQUATION (8)
C...  DOES NOT HAVE ANY DERIVATIVES WITH RESPECT TO Z)
C...
C...          G = HA*F/(HD + HA*F)                                  (13)
C...
C...  EQUATION (13) CAN THEN BE SUBSTITUTED IN EQUATION (7) (WITH
C...  F  = 0)
C...   T
C...
C...          F   = -(HA*F + HD)*HA*F/(HD + HA*F) + HA*F            (14)
C...           ZZ
C...
C...  WHERE THE DERIVATIVE F   IN EQUATION (14) IS NOW AN ORDINARY OR
C...                        ZZ
C...  TOTAL DERIVATIVE SINCE TIME HAS BEEN REMOVED FROM THE PROBLEM.
C...  THUS EQUATION (14) IS A BOUNDARY-VALUE ORDINARY DIFFERENTIAL
C...  EQUATION WHICH IS TO BE INTEGRATED SUBJECT TO BOUNDARY CONDITIONS
C...  (9) AND (10).  TYPICALLY THIS WOULD BE A RELATIVELY DIFFICULT
C...  BOUNDARY-VALUE PROBLEM, BUT NOTE THAT F(Z) = 1 IS A SOLUTION IN
C...  THIS CASE IN THE SENSE THAT IT SATISFIES EQUATIONS (9), (10) AND
C...  (14).  IT ALSO MAKES SENSE PHYSICALLY SINCE AT INFINITE TIME, THE
C...  FLUID IN THE PORE ACHIEVES A UNIFORM CONCENTRATION EQUAL TO THE
C...  CONCENTRATION AT Z = 0.  THUS FOR THE STEADY STATE, THE ENTIRE
C...  PROBLEM HAS THE SIMPLE SOLUTION
C...
C...          F(Z) = 1, G(Z) = HA/(HD + HA)                     (15)(16)
C...
C...  FOR THE NUMERICAL SOLUTION TO BE COMPUTED, HA = HD = H = 1, SO THE
C...  STEADY STATE SOLUTION REDUCES TO
C...
C...          F(Z) = 1, G(Z) = 0.5                              (17)(18)
C...
C...  THE NUMERICAL SOLUTION SHOULD GIVE AN INDICATION OF AGREEMENT
C...  WITH EQUATIONS (17) AND (18) FOR LARGE T.
C...
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
C...
C...  SET THE NUMBER OF GRID POINTS IN Z, INDEX FOR THE PLOTTED SOLUTION
      N=11
      IP=0
C...
C...  SET THE CONSTANTS IN EQUATIONS (7) AND (8)
      H =1.
      HA=1.
      HD=1.
C...
C...  INITIAL CONDITIONS (11) AND (12)
      DO 1 I=1,N
      F(I)=0.
      G(I)=0.
1     CONTINUE
C...
C...  CALL SUBROUTINE DERV TO INITIALZE ALL DERIVATIVES FOR PLOTTING
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
      IF(NCASE.EQ.1)CALL DERV1
      IF(NCASE.EQ.2)CALL DERV2
      IF(NCASE.EQ.3)CALL DERV3
      RETURN
      END
      SUBROUTINE DERV1
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
C...
C...  BOUNDARY CONDITION (9)
       F(1)=1.
      FT(1)=0.
C...
C...  THE PARTIAL DERIVATIVE F   IN EQUATION (7) IS CALCULATED EX-
C...                          ZZ
C...  PLICITY BY SECOND-ORDER CENTRAL DIFFERENCE APPROXIMATIONS
C...
C...  COMPUTE THE SPATIAL INCREMENT IN Z
      DZ=1./FLOAT(N-1)
C...
C...  STEP THROUGH THE N GRID POINTS
      DO 1 I=1,N
C...
C...  AT Z = 0, THE DERIVATIVE NEED NOT BE CALCULATED BECAUSE OF
C...  BOUNDARY CONDITION (9)
      IF(I.EQ.1)THEN
C...
C...  AT Z = 1, THE DERIVATIVE MUST BE CALCULATED TO REFLECT BOUNDARY
C...  CONDITION (10).  IN THIS CASE, A SECOND-ORDER DIFFERENCE APPROXI-
C...  MATION IS APPLIED TO BOUNDARY CONDITION (10) WITH THE RESULT
C...  F(N+1) = F(N-1).  THIS RESULT CAN THEN BE USED TO ELIMINATE F(N+1)
C...  FROM THE FINITE DIFFERENCE APPROXIMATION FOR F   AT Z = 1
C...                                                ZZ
      ELSE IF(I.EQ.N)THEN
      FZZ(N)=2.*(F(N-1)-F(N))/DZ**2
C...
C...  ALL OF THE INTERIOR GRID POINTS CAN BE HANDLED BY THE SAME CENTRAL
C...  DIFFERENCE APPROXIMATION
      ELSE
      FZZ(I)=(F(I+1)-2.*F(I)+F(I-1))/DZ**2
      END IF
1     CONTINUE
C...
C...  PDES
      DO 2 I=1,N
C...
C...     EQUATION (7)
         IF(I.EQ.1)GO TO 3
         FT(I)=FZZ(I)+HD*G(I)-HA*(1.-G(I))*F(I)
C...
C...     EQUATION (8)
3        GT(I)=H*(HA*(1.-G(I))*F(I)-HD*G(I))
2     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
C...
C...  BOUNDARY CONDITION (9)
          F(1)=1.
         FT(1)=0.
C...
C...  F  
C...   Z
      CALL DSS002(0.,1.,N,F,FZ)
C...
C...  BOUNDARY CONDITION (10)
      FZ(N)=0.
C...
C...  F   
C...   ZZ
      CALL DSS002(0.,1.,N,FZ,FZZ)
C...
C...  PDES
      DO 1 I=1,N
C...
C...     EQUATION (7)
         IF(I.EQ.1)GO TO 2
         FT(I)=FZZ(I)+HD*G(I)-HA*(1.-G(I))*F(I)
C...
C...     EQUATION (8)
2        GT(I)=H*(HA*(1.-G(I))*F(I)-HD*G(I))
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
C...
C...  BOUNDARY CONDITION (9)
          F(1)=1.
         FT(1)=0.
C...
C...  F  
C...   Z
      CALL DSS004(0.,1.,N,F,FZ)
C...
C...  BOUNDARY CONDITION (10)
      FZ(N)=0.
C...
C...  F   
C...   ZZ
      CALL DSS004(0.,1.,N,FZ,FZZ)
C...
C...  PDES
      DO 1 I=1,N
C...
C...     EQUATION (7)
         IF(I.EQ.1)GO TO 2
         FT(I)=FZZ(I)+HD*G(I)-HA*(1.-G(I))*F(I)
C...
C...     EQUATION (8)
2        GT(I)=H*(HA*(1.-G(I))*F(I)-HD*G(I))
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/      T,  NSTOP,  NCASE
     1      /Y/  F(11),  G(11)
     2      /F/ FT(11), GT(11)
     3     /SD/ FZ(11),FZZ(11)
     4      /C/      H,     HA,     HD,      N,     IP
C...
C...  DIMENSION THE ARRAYS TO STORE THE NUMERICAL SOLUTION FOR PLOTTING
      DIMENSION   FP(7,11),  GP(7,11),
     1           FTP(7,11), GTP(7,11),
     2           FZP(7,11),FZZP(7,11),
     3            ZP(11)
C...
C...  SELECT THE VALUES OF T AT WHICH THE NUMERICAL SOLUTION IS STORED
C...  FOR SUBSEQUENT PLOTTING
      IF( T.LT.0.005)                  GO TO 1
      IF((T.GT.0.195).AND.(T.LT.0.205))GO TO 1
      IF((T.GT.0.395).AND.(T.LT.0.405))GO TO 1
      IF((T.GT.0.595).AND.(T.LT.0.605))GO TO 1
      IF((T.GT.0.795).AND.(T.LT.0.805))GO TO 1
      IF((T.GT.0.995).AND.(T.LT.1.005))GO TO 1
      IF((T.GT.9.995).AND.(T.LT.10.05))GO TO 1
      RETURN
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING
1     IP=IP+1
      DO 2 I=1,N
        FP(IP,I)=  F(I)
        GP(IP,I)=  G(I)
       FTP(IP,I)= FT(I)
       GTP(IP,I)= GT(I)
       FZP(IP,I)= FZ(I)
      FZZP(IP,I)=FZZ(I)
        ZP(I)=(1./FLOAT(N-1))*FLOAT(I-1)
2     CONTINUE
C...
C...  PRINT THE NUMERICAL SOLUTION AT THE SELECTED VALUES OF T SO THAT
C...  A COMPARISON CAN BE MADE BETWEEN THE THREE CASES
      WRITE(NO,3)T,(ZP(I),I=1,N,2),(F(I),I=1,N,2),(G(I),I=1,N,2)
3     FORMAT(1H ,//,5H T = ,F7.2,/,7H      Z,6F10.3,/,
     1                             7H F(Z,T),6F10.5,/,
     2                             7H G(Z,T),6F10.5)
C...
C...  PLOT THE NUMERICAL SOLUTION AT THE FINAL VALUE OF T
      IF(IP.LT.7)RETURN
      CALL SPLOTS(IP,N,ZP,  FP)
      WRITE(NO,10)
10    FORMAT(///,76H   F(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
      CALL SPLOTS(IP,N,ZP,  GP)
      WRITE(NO,11)
11    FORMAT(///,76H   G(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
      CALL SPLOTS(IP,N,ZP, FTP)
      WRITE(NO,12)
12    FORMAT(///,76H  FT(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
      CALL SPLOTS(IP,N,ZP, GTP)
      WRITE(NO,13)
13    FORMAT(///,76H  GT(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
C...
C...  SINCE THE DERIVATIVES F  AND F   ARE NOT COMPUTED THROUGHOUT THE
C...                         Z      ZZ
C...  DOMAIN 0 LE Z LE 1 IN SUBROUTINE DERV1, THESE DERIVATIVES ARE NOT
C...  PLOTTED IN THE FIRST RUN
      IF(NCASE.EQ.1)RETURN
      CALL SPLOTS(IP,N,ZP, FZP)
      WRITE(NO,14)
14    FORMAT(///,76H  FZ(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
      CALL SPLOTS(IP,N,ZP,FZZP)
      WRITE(NO,15)
15    FORMAT(///,76H FZZ(Z,T) VS Z FOR T = 0.2 (1), 0.4 (2), 0.6 (3), 0.
     18 (4), 1.0 (5), 10.0 (6))
      RETURN
      END
