      SUBROUTINE INITAL
C...
C...  DYNAMIC TWO-DIMENSIONAL TUBULAR REACTOR MODEL IN CYLINDRICAL
C...  COORDINATES
C...
C...  THE MATHEMATICAL MODEL FOR THE REACTOR CONSISTS OF THE FOLLOWING
C...  TWO-DIMENSIONAL PARTIAL DIFFERENTIAL EQUATIONS (PDES)
C...
C...     C  = -V(R)*C  + DC*(C   + (1/R)*C  + C  ) - RK*C        (1)
C...      T          Z        RR          R    ZZ
C...
C...     W  = -V(R)*W  + DW*(W   + (1/R)*W  + W  ) + 2*RK*C      (2)
C...      T          Z        RR          R    ZZ
C...
C...  WITH THE LAMINAR VELOCITY PROFILE
C...
C...     V(R) = 2*VAVG*(1 - R**2/RO**2)
C...
C...  THE CHEMICAL REACTION
C...
C...              KR
C...     E + A  ------+ 2A                                       (3)
C...          2
C...
C...  TAKES PLACE WHERE
C...
C...     E         ELECTRIC DISCHARGE
C...
C...     A         REACTANT
C...      2
C...
C...     A         PRODUCT
C...
C...     KR        FIRST-ORDER REACTION RATE CONSTANT
C...
C...  THE VARIABLES AND PARAMETERS IN EQUATIONS (1) AND (2) ARE
C...
C...     C(R,Z,T)  CONCENTRATION OF A
C...                                 2
C...
C...     W(R,Z,T)  CONCENTRATION OF A
C...
C...     R         RADIAL POSITION WITHIN THE REACTOR
C...
C...     Z         AXIAL POSITION WITHIN THE REACTOR
C...
C...     T         TIME
C...
C...     DC, DW    DIFFUSIVITIES FOR A  AND A, RESPECTIVELY
C...                                  2
C...
C...     VAVG      AVERAGE VELOCITY FOR LAMINAR FLOW IN THE Z
C...               DIRECTION
C...
C...     RO        RADIUS OF THE REACTOR
C...
C...     RK        VOLUMETRIC RATE CONSTANT FOR REACTION (3)
C...
C...  THE INITIAL CONDITIONS FOR EQUATIONS (1) AND (2) ARE
C...
C...     C(R,Z,0) = CI, W(R,Z,0) = 0                          (4)(5)
C...
C...  WHERE
C...
C...     CI        INITIAL CONCENTRATION OF A
C...                                         2
C...
C...  THE BOUNDARY CONDITIONS FOR EQUATION (1) ARE
C...
C...     C (0,Z,T) = C (RO,Z,T) = 0                           (6)(7)
C...      R           R
C...
C...     C(R,0,T) = CI, C (R,ZL,T) = 0                        (8)(9)
C...                     Z
C...
C...  WHERE
C...
C...     CI        ENTERING CONCENTRATION OF A
C...                                          2
C...
C...     ZL        REACTOR LENGTH
C...
C...  THE BOUNDARY CONDITIONS FOR EQUATION (2) ARE
C...
C...     W (0,Z,T) = 0, -DW*W (RO,Z,T) = RD*W(RO,Z,T)**ND   (10)(11)
C...      R                  R
C...
C...     W(R,0,T) = 0, W (R,ZL,T) = 0                       (12)(13)
C...                    Z
C...
C...  WHERE
C...
C...     RD        RATE CONSTANT FOR THE DEPOSITION OF A ON THE
C...               REACTOR WALL
C...
C...     ND        POWER IN THE EQUATION FOR THE RATE OF DESPOSITION
C...               OF A ON THE REACTOR WALL
C...
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/   C(5,25),   W(5,25)
     2      /F/  CT(5,25),  WT(5,25)
     3      /S/  CR(5,25), CRR(5,25),  CZ(5,25), CZZ(5,25),
     4           WR(5,25), WRR(5,25),  WZ(5,25), WZZ(5,25)
     5      /C/        NR,        NZ,        RI,        RO,
     6                 ZO,        ZL,        DR,        DZ,
     7                 DC,        DW,        RK,        RD,
     8               VAVG,        CI,      R(5),      V(5),
     9                 NP,        IP,     Z(25)
      COMMON/IO/       NI,        NO
C...
C...  SET THE NUMBER OF GRID POINTS IN THE RADIAL AND AXIAL
C...  DIRECTIONS
      NR=5
      NZ=25
C...
C...  SET THE REACTOR RADIAL DIMENSIONS
      RI=0.
      RO=1.0
C...
C...  SET THE REACTOR AXIAL DIMENSIONS
      ZO=0.
      ZL=48.
C...
C...  SET THE AVERAGE VELOCITY FOR LAMINAR FLOW, DIFFUSIVITIES
      VAVG=5.
      DC=10.
      DW=14.
C...
C...  SET THE REACTION, DEPOSITION RATE CONSTANTS, POWER FOR THE
C...  ORDER OF DEPOSITION
      RK=0.8
      RD=0.5
      NP=1
C...
C...  COMPUTE THE RADIAL POSITIONS, VELOCITY PROFILE
      DO 1 I=1,NR
      R(I)=(RO/FLOAT(NR-1))*FLOAT(I-1)
      V(I)=2.*VAVG*(1.-R(I)**2/RO**2)
1     CONTINUE
C...
C...  SET THE MODEL INITIAL CONDITIONS
      CI=1.0
      DO 2 J=1,NZ
      Z(J)=ZL*FLOAT(J-1)/(NZ-1)
      DO 2 I=1,NR
      C(I,J)=CI
      W(I,J)=0.
2     CONTINUE
C...
C...  CALCULATE AND PRINT THE INITIAL DERIVATIVES
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/   C(5,25),   W(5,25)
     2      /F/  CT(5,25),  WT(5,25)
     3      /S/  CR(5,25), CRR(5,25),  CZ(5,25), CZZ(5,25),
     4           WR(5,25), WRR(5,25),  WZ(5,25), WZZ(5,25)
     5      /C/        NR,        NZ,        RI,        RO,
     6                 ZO,        ZL,        DR,        DZ,
     7                 DC,        DW,        RK,        RD,
     8               VAVG,        CI,      R(5),      V(5),
     9                 NP,        IP,     Z(25)
C...
C...  CONSTRAIN THE SOLUTION TO MEANINGFUL VALUES
      DO 8 I=1,NR
      DO 8 J=1,NZ
      IF(C(I,J).LT.0.)C(I,J)=0.
      IF(W(I,J).LT.0.)W(I,J)=0.
8     CONTINUE
C...
C...  SET THE BOUNDARY CONDITIONS AT Z = 0
      DO 1 I=1,NR
      C(I,1)=CI
      W(I,1)=0.
      CT(I,1)=0.
      WT(I,1)=0.
1     CONTINUE
C...
C...  COMPUTE THE FIRST-ORDER DERIVATIVES IN R
      CALL DSS034(RI,RO,NR,NZ,1,C,CR,0.)
      CALL DSS034(RI,RO,NR,NZ,1,W,WR,0.)
C...
C...  SET THE BOUNDARY CONDITIONS AT R = RI
      DO 2 J=1,NZ
      CR(1,J)=0.
      WR(1,J)=0.
2     CONTINUE
C...
C...  SET THE BOUNDARY CONDITIONS AT R = RO
      DO 3 J=1,NZ
      CR(NR,J)=0.
      WR(NR,J)=0.
      IF(W(NR,J).GT.0.)WR(NR,J)=-(1./DW)*W(NR,J)**NP
3     CONTINUE
C...
C...  COMPUTE THE SECOND-ORDER DERIVATIVES IN R
      CALL DSS034(RI,RO,NR,NZ,1,CR,CRR,0.)
      CALL DSS034(RI,RO,NR,NZ,1,WR,WRR,0.)
C...
C...  COMPUTE THE FIRST-ORDER DERIVATIVES IN Z FOR DIFFUSION
      CALL DSS034(ZO,ZL,NR,NZ,2,C,CZ,0.)
      CALL DSS034(ZO,ZL,NR,NZ,2,W,WZ,0.)
C...
C...  SET THE BOUNDARY CONDITIONS AT Z = ZL
      DO 4 I=1,NR
      CZ(I,NZ)=0.
      WZ(I,NZ)=0.
4     CONTINUE
C...
C...  COMPUTE THE SECOND-ORDER DERIVATIVES IN Z FOR DIFFUSION
      CALL DSS034(ZO,ZL,NR,NZ,2,CZ,CZZ,0.)
      CALL DSS034(ZO,ZL,NR,NZ,2,WZ,WZZ,0.)
C...
C...  COMPUTE THE FIRST-ORDER DERIVATIVES IN Z FOR CONVECTION
      CALL DSS034(ZO,ZL,NR,NZ,2,C,CZ,1.)
      CALL DSS034(ZO,ZL,NR,NZ,2,W,WZ,1.)
C...
C...  SET THE BOUNDARY CONDITIONS AT Z = ZL
      DO 5 I=1,NR
      CZ(I,NZ)=0.
      WZ(I,NZ)=0.
5     CONTINUE
C...
C...  ASSEMBLE EQUATION (1)
C...  ******************************************************************
      DO 6 I=1,NR
      DO 6 J=2,NZ
      RATE=0.
C...
C...  INCLUDE THE REACTION ZONE
      IF(J.EQ.11)RATE=-RK*C(I,J)/2.
      IF(J.EQ.16)RATE=-RK*C(I,J)/2.
      IF((J.GT.11).AND.(J.LT.16))RATE=-RK*C(I,J)
C...
C...  INDETERMINATE FORM AT R = 0
      IF(I.EQ.1)
     1 CT(I,J)=-V(I)*CZ(I,J)+DC*(2.*CRR(I,J)+CZZ(I,J))+RATE
C...
C...  INCLUDE THE TERM (1/R)*C  FOR R NE 0
C...                          R
      IF(I.GT.1)
     1 CT(I,J)=-V(I)*CZ(I,J)+DC*(CRR(I,J)+(1./R(I))*CR(I,J)+CZZ(I,J))
     2         +RATE
6     CONTINUE
C...
C...  ASSEMBLE EQUATION (2)
C...  ******************************************************************
      DO 7 I=1,NR
      DO 7 J=2,NZ
      RATE=0.
C...
C...  INCLUDE THE REACTION ZONE
      IF(J.EQ.11)RATE=RK*C(I,J)
      IF(J.EQ.16)RATE=RK*C(I,J)
      IF((J.GT.11).AND.(J.LT.16))RATE=2.*RK*C(I,J)
C...
C...  INDETERMINATE FORM AT R = 0
      IF(I.EQ.1)
     1 WT(I,J)=-V(I)*WZ(I,J)+DW*(2.*WRR(I,J)+WZZ(I,J))+RATE
C...
C...  INCLUDE THE TERM (1/R)*W  FOR R NE 0
C...                          R
      IF(I.GT.1)
     1 WT(I,J)=-V(I)*WZ(I,J)+DW*(WRR(I,J)+(1./R(I))*WR(I,J)+WZZ(I,J))
     2         +RATE
7     CONTINUE
C...  ******************************************************************
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,      NFIN,     NORUN
     1      /Y/   C(5,25),   W(5,25)
     2      /F/  CT(5,25),  WT(5,25)
     3      /S/  CR(5,25), CRR(5,25),  CZ(5,25), CZZ(5,25),
     4           WR(5,25), WRR(5,25),  WZ(5,25), WZZ(5,25)
     5      /C/        NR,        NZ,        RI,        RO,
     6                 ZO,        ZL,        DR,        DZ,
     7                 DC,        DW,        RK,        RD,
     8               VAVG,        CI,      R(5),      V(5),
     9                 NP,        IP,     Z(25)
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,1)T,((C(I,J),J=1,NZ,4),I=1,NR)
1     FORMAT(1H ,//,5H t = ,F7.2,/,9H C(r,z,t),/,
     + 10X,8H  z = 0 ,8H  z = 8 ,8H z = 16 ,
     +     8H z = 24 ,8H z = 32 ,8H z = 40 ,
     +     8H z = 48 ,/,
     +  9H r = 1   ,7F8.4,/,
     +  9H r = 0.75,7F8.4,/,
     +  9H r = 0.5 ,7F8.4,/,
     +  9H r = 0.25,7F8.4,/,
     +  9H r = 0   ,7F8.4,/)
      WRITE(NO,2)  ((W(I,J),J=1,NZ,4),I=1,NR)
2     FORMAT(1H ,/,                9H W(r,z,t),/,
     + 10X,8H  z = 0 ,8H  z = 8 ,8H z = 16 ,
     +     8H z = 24 ,8H z = 32 ,8H z = 40 ,
     +     8H z = 48 ,/,
     +  9H r = 1   ,7F8.4,/,
     +  9H r = 0.75,7F8.4,/,
     +  9H r = 0.5 ,7F8.4,/,
     +  9H r = 0.25,7F8.4,/,
     +  9H r = 0   ,7F8.4,/)
C...
C...  PRINT THE INITIAL DERIVATIVES
      IF(IP.NE.-1)GO TO 100
      WRITE(NO,3)(( CR(I,J),J=1,NZ,4),I=1,NR)
3     FORMAT(1H , /,           11H  CR(R,Z,T),/,(7F8.4))
      WRITE(NO,4)((CRR(I,J),J=1,NZ,4),I=1,NR)
4     FORMAT(1H , /,           11H CRR(R,Z,T),/,(7F8.4))
      WRITE(NO,5)(( CZ(I,J),J=1,NZ,4),I=1,NR)
5     FORMAT(1H , /,           11H  CZ(R,Z,T),/,(7F8.4))
      WRITE(NO,6)((CZZ(I,J),J=1,NZ,4),I=1,NR)
6     FORMAT(1H , /,           11H CZZ(R,Z,T),/,(7F8.4))
      WRITE(NO,7)(( CT(I,J),J=1,NZ,4),I=1,NR)
7     FORMAT(1H , /,           11H  CT(R,Z,T),/,(7F8.4))
      WRITE(NO,8)(( WR(I,J),J=1,NZ,4),I=1,NR)
8     FORMAT(1H , /,           11H  WR(R,Z,T),/,(7F8.4))
      WRITE(NO,9)((WRR(I,J),J=1,NZ,4),I=1,NR)
9     FORMAT(1H , /,           11H WRR(R,Z,T),/,(7F8.4))
      WRITE(NO,10)(( WZ(I,J),J=1,NZ,4),I=1,NR)
10    FORMAT(1H , /,           11H  WZ(R,Z,T),/,(7F8.4))
      WRITE(NO,11)((WZZ(I,J),J=1,NZ,4),I=1,NR)
11    FORMAT(1H , /,           11H WZZ(R,Z,T),/,(7F8.4))
      WRITE(NO,12)(( WT(I,J),J=1,NZ,4),I=1,NR)
12    FORMAT(1H , /,           11H  WT(R,Z,T),/,(7F8.4))
100   IP=IP+1
      RETURN
      END
