      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A PACKED, ADIABATIC TUBULAR REACTOR
C...
C...  THE FOLLOWING MODEL OF THE DYNAMICS OF A PACKED, ADIABATIC
C...  TUBULAR REACTOR IS PRESENTED BY LEE, E. STANLEY, QUASILINEAR-
C...  IZATION AND INVARIANT IMBEDDING, ACADEMIC PRESS, NY, 1968, PP
C...  285-292
C...
C...     P  = -P  + (1/PE)*P   - R*(P**2)*EXP(-E/(R*TK))         (1)
C...      T     Z           ZZ
C...
C...     TK  = -TK  + (1/PE)*TK   - Q*R*(P**2)*EXP(-E/(R*TK))    (2)
C...       T      Z            ZZ
C...
C...  BWHERE
C...
C...     P        PARTIAL PRESSURE OF THE REACTANT (ATM)
C...
C...     TK       TEMPERATURE OF FLUID STREAM (K)
C...
C...     T        TIME (DIMENSIONLESS)
C...
C...     Z        POSITION ALONG THE REACTOR (DIMENSIONLESS)
C...
C...     PE       PECLET NUMBER (DP*V/D)
C...
C...     DP       PACKING PARTICLE DIAMETER (CM)
C...
C...     D        AXIAL DISPERSION COEFFICIENT
C...
C...     V        LINEAR VELOCITY (CM/SEC)
C...
C...     R        REACTION GROUP (K*DP/V)
C...
C...     Q        REACTION GROUP (DH*RHOM/(PT*MW*RHOF*CPF))
C...
C...     K        REACTION SPECIFIC RATE CONSTANT (1/SEC-ATM)
C...
C...     DH       HEAT OF REACTION (CAL/GM)
C...
C...     RHOM     FLUID MOLAR DENSITY (GM MOL/CM**3)
C...
C...     PT       TOTAL PRESSURE (ATM)
C...
C...     MW       MOLECULAR WEIGHT OF THE REACTANT (GM/GM MOL)
C...
C...     RHOF     FLUID DENSITY (GM/CM**3)
C...
C...     CPF      FLUID SPECIFIC HEAT (CAL/GM-K)
C...
C...     E        ACTIVATION ENERGY (CAL/GM MOL)
C...
C...     R        GAS CONSTANT (CAL/GM MOL-K)
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATION (1) ARE
C...
C...     P(Z,0) = 0                                              (3)
C...
C...     PEN = P(0,T) - (1/PE)*P (0,T), P (ZL,T) = 0          (4)(5)
C...                            Z        Z
C...  WHERE
C...
C...     PEN      ENTERING REACTANT PARTIAL PRESSURE (ATM)
C...
C...     ZL       REACTOR LENGTH (DIMENSIONLESS)
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATION (2) ARE
C...
C...     TK(Z,0) = 0                                             (6)
C...
C...     TEN = TK(0,T) - (1/PE)*TK (0,T), TK (ZL,T) = 0       (7)(8)
C...                              Z         Z
C...  WHERE
C...
C...     TEN      ENTERING REACTANT TEMPERATURE (K)
C...
C...  NOTE THAT Z AND T EQUAL DIMENSIONAL TIMES NORMALZED BY DP AND 
C...  DP**2/D, RESPECTIVELY, I.E, Z = Z/DP, T = T*D/DP**2 WHERE THE
C...  LHS Z AND T ARE DIMENSIONLESS.
C...
C...  LEE REPORTS THE FOLLOWING NUMERICAL RESULT FROM A QUASILINEARI-
C...  ZATION/FINITE DIFFERENCE SOLUTION (TABLES 9.2A, 9.2B, P291): 
C...  P(0,1) = 0.057121, TK(0,1) = 1254.1 (OR 1254.1/1.8 = 696.7), 
C...  FOR THE PARAMETERS PE = 2, PEN = 0.07, TEN = 694.4, R = 0.5E+08, 
C...  Q = -555.5, E/R = 12,222, ZF = 48.
C...
C...  THE FOLLOWING PROGRAM MAY NOT REPRODUCE THIS RESULT PRECISELY
C...  FOR AT LEAST TWO REASONS:
C...  
C...     (1)  LEE USED A SPATIAL INCREMENT OF 0.1 WHILE THE FOLLOWING
C...          PROGRAM USES 1 (NZ = 49, ZF = 48, DZ = 48/(49-1) = 1).
C...
C...     (2)  LEE USED SECOND ORDER APPROXIMATIONS FOR THE SPATIAL
C...          DERIVATIVES, WHILE SUBROUTINES DSS020 AND DSS044 ARE
C...          BASED ON FOURTH ORDER DIFFERENCES (WHICH WILL TEND TO
C...          OFFSET THE SMALLER SPATIAL INTERVAL USED BY LEE).
C...
C...  TO INVESTIGATE POINT (2), IN PARTICULAR, THE SPATIAL CONVERGENCE
C...  OF THE SOLUTION, THE FOLLOWING PROGRAM COMPUTES SOLUTIONS FOR NZ
C...  = 49 AND 98.  THE OUTPUT OF THE PROGRAM FOR Z = 0, T = 1 IS
C
C                    NZ=49                        NZ=96
C
C           P(Z=0,T=1)   TK(Z=0,T=1)     P(Z=0,T=1)   TK(Z=0,T=1)
C
C             0.05775       696.6          0.05875       696.5
C
C     THE FOLLOWING IS A SUMMARY OF THE OUTPUT REPORTED BY LEE (P292) FOR
C     Z = 48, T = 60 AND 70
C
C      P(Z=48,T=60)   TK(Z=48,T=60)     P(Z=48,T=70)   TK(Z=48,T=70)
C
C        0.0088082     1308.9/1.8 =       0.0086671     1311.2/1.8 =
C                         727.2                            728.4
C
C     THE CORRESPONDING OUTPUT FROM THE FOLLOWING PROGRAM IS:
C
C           P(Z=48,T=60)   TK(Z=48,T=60)     P(Z=48,T=70)   TK(Z=48,T=70)
C
C     NZ=49   0.00886          726.9           0.00870          728.3 
C
C     NZ=98   0.00887          726.9           0.00871          728.3 
C
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  NORUN = 1, NZ = 49
      IF(NORUN.EQ.1)CALL INIT49
C...
C...  NORUN = 1, NZ = 98
      IF(NORUN.EQ.2)CALL INIT98
      RETURN
      END
      SUBROUTINE INIT49
      PARAMETER (NZ=49)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  MODEL PARAMETERS
      PE=2.0
      PEN=0.07
      TEN=694.4
      R=0.5E+08
      ZF=48.0
      Q=-555.5
      ER=12222.0
C...
C...  SPATIAL GRID
      DO 1 I=1,NZ
         ZG(I)=ZF*FLOAT(I-1)/FLOAT(NZ-1)
1     CONTINUE
C...
C...  INITIAL CONDITIONS
      DO 2 I=1,NZ
         P(I)=0.
        TK(I)=705.5
2     CONTINUE
C...
C...  INITIAL DERIVATIVE
      CALL DERV
C...
C...  INITIALIZE SOME COUNTERS USED IN SUBROUTINE PRINT
      IP=0
      IS=0
      RETURN
      END      
      SUBROUTINE INIT98
      PARAMETER (NZ=98)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  MODEL PARAMETERS
      PE=2.0
      PEN=0.07
      TEN=694.4
      R=0.5E+08
      ZF=48.0
      Q=-555.5
      ER=12222.0
C...
C...  SPATIAL GRID
      DO 1 I=1,NZ
         ZG(I)=ZF*FLOAT(I-1)/FLOAT(NZ-1)
1     CONTINUE
C...
C...  INITIAL CONDITIONS
      DO 2 I=1,NZ
         P(I)=0.
        TK(I)=705.5
2     CONTINUE
C...
C...  INITIAL DERIVATIVE
      CALL DERV
C...
C...  INITIALIZE SOME COUNTERS USED IN SUBROUTINE PRINT
      IP=0
      IS=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  NORUN = 1, NZ = 49
      IF(NORUN.EQ.1)CALL DERV49
C...
C...  NORUN = 1, NZ = 98
      IF(NORUN.EQ.2)CALL DERV98
      RETURN
      END
      SUBROUTINE DERV49
      PARAMETER (NZ=49)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  P
C...   Z
      CALL DSS020(0.,ZF,NZ,P,PZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=2
      PZ(1)=PE*(P(1)-PEN)
C...
C...  BOUNDARY CONDITION AT Z = ZF
      NU=2
      PZ(NZ)=0.
C...
C...  P
C...   ZZ
      CALL DSS044(0.,ZF,NZ,P,PZ,PZZ,NL,NU)
C...
C...  TK
C...    Z
      CALL DSS020(0.,ZF,NZ,TK,TKZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=2
      TKZ(1)=PE*(TK(1)-TEN)
C...
C...  BOUNDARY CONDITION AT Z = ZF
      NU=2
      TKZ(NZ)=0.
C...
C...  TK
C...    ZZ
      CALL DSS044(0.,ZF,NZ,TK,TKZ,TKZZ,NL,NU)
C...
C...  PDES
      DO 1 I=1,NZ
         REACT=R*(P(I)**2)*EXP(-ER/TK(I))
          PT(I)=- PZ(I)+(1.0/PE)* PZZ(I)-  REACT
         TKT(I)=-TKZ(I)+(1.0/PE)*TKZZ(I)-Q*REACT
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV98
      PARAMETER (NZ=98)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  P
C...   Z
      CALL DSS020(0.,ZF,NZ,P,PZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=2
      PZ(1)=PE*(P(1)-PEN)
C...
C...  BOUNDARY CONDITION AT Z = ZF
      NU=2
      PZ(NZ)=0.
C...
C...  P
C...   ZZ
      CALL DSS044(0.,ZF,NZ,P,PZ,PZZ,NL,NU)
C...
C...  TK
C...    Z
      CALL DSS020(0.,ZF,NZ,TK,TKZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=2
      TKZ(1)=PE*(TK(1)-TEN)
C...
C...  BOUNDARY CONDITION AT Z = ZF
      NU=2
      TKZ(NZ)=0.
C...
C...  TK
C...    ZZ
      CALL DSS044(0.,ZF,NZ,TK,TKZ,TKZZ,NL,NU)
C...
C...  PDES
      DO 1 I=1,NZ
         REACT=R*(P(I)**2)*EXP(-ER/TK(I))
          PT(I)=- PZ(I)+(1.0/PE)* PZZ(I)-  REACT
         TKT(I)=-TKZ(I)+(1.0/PE)*TKZZ(I)-Q*REACT
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  NORUN = 1, NZ = 49
      IF(NORUN.EQ.1)CALL PRNT49(NI,NO)
C...
C...  NORUN = 1, NZ = 98
      IF(NORUN.EQ.2)CALL PRNT98(NI,NO)
      RETURN
      END
      SUBROUTINE PRNT49(NI,NO)
      PARAMETER (NZ=49)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  ARRAY FOR THE PLOTTED SOLUTION
      PARAMETER (NS=8)
      REAL PP(NS,NZ), TKP(NS,NZ)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)THEN
         WRITE(NO,1)
1        FORMAT(9X,'T',5X,'ZG(1)',5X,' P(1)',4X,'ZG(NZ)',4X,' P(NZ)',/,
     +             10X,5X,'ZG(1)',5X,'TK(1)',4X,'ZG(NZ)',4X,'TK(NZ)')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,ZG(1), P(1),ZG(NZ), P(NZ),
     +             ZG(1),TK(1),ZG(NZ),TK(NZ)
2     FORMAT(F10.0,F10.0,F10.5,F10.0,F10.5,/,
     +         10X,F10.0,F10.1,F10.0,F10.1,/)
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING AT T = 5, 10, 20, 30,
C...  40, 50, 60, 70
      IF((IP.EQ. 5).OR.(IP.EQ.10).OR.(IP.EQ.20).OR.
     +   (IP.EQ.30).OR.(IP.EQ.40).OR.(IP.EQ.50).OR.
     +   (IP.EQ.60).OR.(IP.EQ.70))THEN
         IS=IS+1
         DO 3 I=1,NZ
            PP(IS,I)=P(I)
           TKP(IS,I)=TK(I)
3        CONTINUE
      END IF
      IP=IP+1
C...
C...  AT THE END OF THE SOLUTION, T = 70, PLOT THE SPATIAL PROFILES
      IF(IP.LT.71)RETURN
C...
C...  PARTIAL PRESSURE PROFILES
      CALL TPLOTS(IS,NZ,ZG,PP)
      WRITE(NO,4)
4     FORMAT(/,' P(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60, 70')
C...
C...  TEMPERATURE PROFILES
      CALL TPLOTS(IS,NZ,ZG,TKP)
      WRITE(NO,5)
5     FORMAT(/,' TK(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60, 70')
      RETURN
      END
      SUBROUTINE PRNT98(NI,NO)
      PARAMETER (NZ=98)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ),    TK(NZ)
     +      /F/    PT(NZ),   TKT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),   TKZ(NZ),  TKZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,       TEN,         R,        ZF,
     +                  Q,        ER
     +      /I/        IP,        IS
C...
C...  ARRAY FOR THE PLOTTED SOLUTION
      PARAMETER (NS=8)
      REAL PP(NS,NZ), TKP(NS,NZ)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)THEN
         WRITE(NO,1)
1        FORMAT(9X,'T',5X,'ZG(1)',5X,' P(1)',4X,'ZG(NZ)',4X,' P(NZ)',/,
     +             10X,5X,'ZG(1)',5X,'TK(1)',4X,'ZG(NZ)',4X,'TK(NZ)')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,ZG(1), P(1),ZG(NZ), P(NZ),
     +             ZG(1),TK(1),ZG(NZ),TK(NZ)
2     FORMAT(F10.0,F10.0,F10.5,F10.0,F10.5,/,
     +         10X,F10.0,F10.1,F10.0,F10.1,/)
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING AT T = 5, 10, 20, 30,
C...  40, 50, 60, 70
      IF((IP.EQ. 5).OR.(IP.EQ.10).OR.(IP.EQ.20).OR.
     +   (IP.EQ.30).OR.(IP.EQ.40).OR.(IP.EQ.50).OR.
     +   (IP.EQ.60).OR.(IP.EQ.70))THEN
         IS=IS+1
         DO 3 I=1,NZ
            PP(IS,I)=P(I)
           TKP(IS,I)=TK(I)
3        CONTINUE
      END IF
      IP=IP+1
C...
C...  AT THE END OF THE SOLUTION, T = 70, PLOT THE SPATIAL PROFILES
      IF(IP.LT.71)RETURN
C...
C...  PARTIAL PRESSURE PROFILES
      CALL TPLOTS(IS,NZ,ZG,PP)
      WRITE(NO,4)
4     FORMAT(/,' P(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60, 70')
C...
C...  TEMPERATURE PROFILES
      CALL TPLOTS(IS,NZ,ZG,TKP)
      WRITE(NO,5)
5     FORMAT(/,' TK(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60, 70')
      RETURN
      END
