      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A PACKED, ISOTHERMAL TUBULAR REACTOR
C...
C...  THE FOLLOWING MODEL OF THE DYNAMICS OF A PACKED, ISOTHERMAL
C...  TUBULAR REACTOR IS PRESENTED BY LEE, E. STANLEY, QUASILINEAR-
C...  IZATION AND INVARIANT IMBEDDING, ACADEMIC PRESS, NY, 1968, PP
C...  279-285
C...
C...     P  = -P  + (1/PE)*P   - R*(P**2)                        (1)
C...      T     Z           ZZ
C...
C...  WHERE
C...
C...     P        PARTIAL PRESSURE OF THE REACTANT (ATM)
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...     K        REACTION SPECIFIC RATE CONSTANT (1/SEC-ATM)
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATION (1) ARE
C...
C...     P(Z,0) = 0                                              (2)
C...
C...     PEN = P(0,T) - (1/PE)*P (0,T), P (ZL,T) = 0          (3)(4)
C...                            Z        Z
C...  WHERE
C...
C...     PEN      ENTERING REACTANT PARTIAL PRESSURE (ATM)
C...
C...     ZL       REACTOR LENGTH (DIMENSIONLESS)
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 (TABLE 9.1, P284): P(0,1) = 
C...  0.057238 FOR THE PARAMETERS PE = 2, PEN = 0.07, R = 1, ZF = 48.
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.
C...
C...  THE FOLLOWING IS A SUMMARY OF THE OUTPUT
C...
C...                     P(Z=0,T=1)     P(Z=0,T=60)    P(Z=48,T=60)
C...
C...     NZ = 49          0.05786         0.06783        0.01624
C...
C...     NZ = 98          0.05889         0.06783        0.01626
C...
C...  THESE RESULTS SUGGEST THE STEADY STATE SOLUTION (T = 60) HAS CONVER-
C...  GED TO ABOUT FOUR FIGURES.
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)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /I/        IP,        IS
C...
C...  MODEL PARAMETERS
      PE=2.0
      PEN=0.07
      R=1.0
      ZF=48.0
C...
C...  SPATIAL GRID
      DO 1 I=1,NZ
         ZG(I)=ZF*FLOAT(I-1)/FLOAT(NZ-1)
1     CONTINUE
C...
C...  INITIAL CONDITION
      DO 2 I=1,NZ
         P(I)=0.
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)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /I/        IP,        IS
C...
C...  MODEL PARAMETERS
      PE=2.0
      PEN=0.07
      R=1.0
      ZF=48.0
C...
C...  SPATIAL GRID
      DO 1 I=1,NZ
         ZG(I)=ZF*FLOAT(I-1)/FLOAT(NZ-1)
1     CONTINUE
C...
C...  INITIAL CONDITION
      DO 2 I=1,NZ
         P(I)=0.
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)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /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...  PDE
      DO 1 I=1,NZ
         PT(I)=-PZ(I)+(1.0/PE)*PZZ(I)-R*P(I)**2
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV98
      PARAMETER (NZ=98)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /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...  PDE
      DO 1 I=1,NZ
         PT(I)=-PZ(I)+(1.0/PE)*PZZ(I)-R*P(I)**2
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)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /I/        IP,        IS
C...
C...  ARRAY FOR THE PLOTTED SOLUTION
      PARAMETER (NS=7)
      REAL PP(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)',6X,'P(1)',4X,'ZG(NZ)',5X,'P(NZ)')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,ZG(1),P(1),ZG(NZ),P(NZ)
2     FORMAT(F10.0,F10.0,F10.5,F10.0,F10.5)
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING AT T = 5, 10, 20, 30,
C...  40, 50, 60
      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))THEN
         IS=IS+1
         DO 3 I=1,NZ
            PP(IS,I)=P(I)
3        CONTINUE
      END IF
      IP=IP+1
C...
C...  AT THE END OF THE SOLUTION, T = 60, PLOT THE SPATIAL PROFILES
      IF(IP.LT.61)RETURN
      CALL TPLOTS(IS,NZ,ZG,PP)
      WRITE(NO,4)
4     FORMAT(/,' P(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60')
      RETURN
      END
      SUBROUTINE PRNT98(NI,NO)
      PARAMETER (NZ=98)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/     P(NZ)
     +      /F/    PT(NZ)
     +      /S/    PZ(NZ),   PZZ(NZ),    ZG(NZ)
     +      /C/        PE,       PEN,         R,        ZF
     +      /I/        IP,        IS
C...
C...  ARRAY FOR THE PLOTTED SOLUTION
      PARAMETER (NS=7)
      REAL PP(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)',6X,'P(1)',4X,'ZG(NZ)',5X,'P(NZ)')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,ZG(1),P(1),ZG(NZ),P(NZ)
2     FORMAT(F10.0,F10.0,F10.5,F10.0,F10.5)
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING AT T = 5, 10, 20, 30,
C...  40, 50, 60
      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))THEN
         IS=IS+1
         DO 3 I=1,NZ
            PP(IS,I)=P(I)
3        CONTINUE
      END IF
      IP=IP+1
C...
C...  AT THE END OF THE SOLUTION, T = 60, PLOT THE SPATIAL PROFILES
      IF(IP.LT.61)RETURN
      CALL TPLOTS(IS,NZ,ZG,PP)
      WRITE(NO,4)
4     FORMAT(/,' P(Z,T) VS Z, T = 5, 10, 20, 30, 40, 50, 60')
      RETURN
      END
