      SUBROUTINE INITAL
C...
C...  NONISOTHERMAL CATALYTIC TUBULAR REACTOR
C...
C...  DEVELOP A MATHEMATICAL MODEL AND SIMULATION TO STUDY THE RESPONSE
C...  OF AN ADIABATIC, GAS-PHASE, PACKED TUBULAR REACTOR.  THE REACTOR
C...  IS A CYLINDRICAL VESSEL, 1 M ID BY 6 M LONG, PACKED WITH A METAL
C...  CATALYST THAT OCCUPIES 5 PERCENT OF THE TOTAL VOLUME, PROVIDES
C...  150 M**2 OF AREA PER M**3 OF TOTAL VOLUME, WEIGHTS 6000 KG/M**3
C...  AND HAS A HEAT CAPACITY OF 400 J/(KG-C).  THE HEAT TRANSFER CO-
C...  EFFICIENT BETWEEN THE CATALYST AND THE GAS IS 60 J/(SEC-M**3-C).
C...  THE PRESSURE AT THE INLET OF THE REACTOR IS 7.0E+05 PA AND THE
C...  PRESSURE DROP OVER THE LENGTH OF THE REACTOR IS 3.0E+04 PA AT
C...  THE DESIGN SUPERFICIAL VELOCITY OF 0.3 M/SEC AT INLET CONDITIONS.
C...
C...  THE REACTION IS FIRST ORDER
C...
C...         K
C...     A -----+ B
C...
C...  THE HEAT OF REACTION IS -70000 J/MOLE AND THE REACTION RATE CON-
C...  STANT IS GIVEN BY K = 12000*EXP(-7500/TC) (1/SEC) (TC IN K).  ASSUME
C...  PERFECT GASES WITH MOLECULAR WEIGHTS OF 40 AND HEAT CAPACITIES OF
C...  OF 600 J/(KG-C).
C...
C...  THE MODEL FOR THE NONISOTHERMAL CATALYTIC TUBULAR REACTOR CODED
C...  IN SUBROUTINES INITAL, DERV AND PRINT INCLUDES
C...
C...  (1)  ENERGY BALANCE FOR THE CATALYST
C...       -------------------------------
C...
C...       RATE OF ACCUMULATION OF ENERGY WITHIN THE CATALYST =
C...
C...       HEAT TRANSFER RATE TO THE CATALYST
C...
C...       (S*DZ*(1-EPS)*RHOC*CPC*TC)  = H*AV*S*DZ*(TF - TC)
C...                                 T
C...
C...       DIVISION BY S*DZ*(1-EPS)*RHOC*CPC GIVES
C...
C...       TC  = H*AV/((1-EPS)*RHOC*CPC)*(TF - TC)                  (1)
C...         T
C...
C...       CHECK ON UNITS
C...
C...       (S*DZ*(1-EPS)*RHOC*CPC*TC)
C...                                 T
C...
C...       (M**2)*(M)*(KG/M**3)*(J/(KG-C))*(C)*(1/SEC) = J/SEC
C...
C...       H*AV*S*DZ*(TF - TC)
C...
C...       (J/(SEC-M**2-C))*(M**2/M**3)*(M**2)*(M)*(C) = J/SEC
C...
C...  (2)  ENERGY BALANCE FOR THE FLUID
C...       ----------------------------
C...
C...       RATE OF ACCUMULATION OF ENERGY WITHIN THE FLUID =
C...
C...       NET FLOW RATE OF ENERGY OF THE FLUID +
C...
C...       HEAT TRANSFER RATE TO THE FLUID +
C...
C...       RATE OF HEAT PRODUCTION BY THE REACTION
C...
C...       (S*DZ*EPS*RHOF*CVF*TF)  =
C...                             T
C...
C...       (V*S*RHOF*CPF*TF)     - (V*S*RHOF*CPF*TF)
C...                        AT Z                    AT Z+DZ
C...
C...       - H*AV*S*DZ*(TF - TC)
C...
C...       - S*DZ*EPS*DH*K*CA
C...
C...       DIVISION BY (S*DZ*EPS*CVF) WITH DZ ----+ 0 GIVES
C...
C...       (RHOF*TF)  = - (CPF/(EPS*CVF))*(RHOF*V*TF)
C...                T                                Z
C...                                                                (2)
C...                    - H*AV/(EPS*CVF)*(TF - TC) - (DH*K*CA)/CVF
C...
C...       CHECK ON UNITS
C...
C...       (S*DZ*EPS*RHOF*CVF*TF)
C...                             T
C...
C...       (M**2)*(M)*(KG/M**3)*(J/(KG-C))*(C)*(1/SEC) = J/SEC
C...
C...       (V*S*RHOF*CPF*TF)
C...
C...       (M/SEC)*(M**2)*(KG/M**3)*(J/(KG-C))*(C) = J/SEC
C...
C...       H*AV*S*DZ*(TF - TC)
C...
C...       (J/(SEC-M**2-C))*(M**2/M**3)*(M**2)*(M)*(C) = J/SEC
C...
C...       S*DZ*EPS*DH*K*CA
C...
C...       (M**2)*(M)*(J/MOLE)*(1/SEC)*(MOLES/M**3) = J/SEC
C...
C...  (3)  TOTAL MASS BALANCE FOR THE FLUID
C...       --------------------------------
C...
C...       RATE OF ACCUMULATION OF MASS WITHIN THE FLUID =
C...
C...       NET FLOW RATE OF MASS WITHIN THE FLUID
C...
C...       (S*DZ*EPS*RHOF)  = (V*S*RHOF)     - (V*S*RHOF)
C...                      T             AT Z             AT Z+DZ
C...
C...       DIVISION BY (S*DZ*EPS) WITH DZ ----+ 0 GIVES
C...
C...       RHOF  = -(1/EPS)*(V*RHOF)                                (3)
C...           T                   Z
C...
C...       CHECK ON UNITS
C...
C...       (S*DZ*EPS*RHOF)
C...                      T
C...
C...       (M**2)*(M)*(KG/MM**3)*(1/SEC) = KG/SEC
C...
C...       V*S*RHOF
C...
C...       (M/SEC)*(M**2)*(KG/M**3) = KG/SEC
C...
C...  (4)  COMPONENT BALANCE FOR THE FLUID
C...       -------------------------------
C...
C...       RATE OF ACCUMULATION OF A WITHIN THE FLUID =
C...
C...       NET FLOW OF A WITHIN THE FLUID -
C...
C...       RATE OF CONSUMPTION OF A WITHIN THE FLUID BY REACTION
C...
C...       (S*DZ*EPS*CA)  =
C...                    T
C...
C...       (V*S*CA)     - (V*S*CA)
C...               AT Z           AT Z+DZ
C...
C...       - S*DZ*EPS*K*CA
C...
C...       DIVISION BY (S*DZ*EPS) WITH DZ ----+ 0 GIVES
C...
C...       CA  = -(1/EPS)*(V*CA)  - K*CA                            (4)
C...         T                  Z
C...
C...       CHECK ON UNITS
C...
C...       (S*DZ*EPS*CA)
C...                    T
C...
C...       (M**2)*(M)*(MOLE/M**3)*(1/SEC) = MOLE/SEC
C...
C...       (V*S*CA)
C...
C...       (M/SEC)*(M**2)*(MOLE/M**3) = MOLE/SEC
C...
C...       (S*DZ*EPS*K*CA)
C...
C...       (M**2)*(M)*(1/SEC)*(MOLE/M**3) = MOLE/SEC
C...
C...  (5)  MOMENTUM BALANCE FOR THE FLUID
C...       ------------------------------
C...
C...       RATE OF ACCUMULATION OF MOMENTUM WITHIN THE FLUID =
C...       (SINCE THIS IS NEWTON*S SECOND LAW, THIS TERM
C...        COULD BE CALLED THE INERTIAL FORCE OR THE TIME
C...        RATE OF CHANGE OF MOMENTUM.  STILL ANOTHER WAY
C...        TO EXPRESS THIS TERM IS TO SAY THAT IT IS THE
C...        M*A IN F = M*A, IN WHICH CASE, THE RHS TERMS OF
C...        THE EQUATION ARE THE NET FORCES ACTING ON THE
C...        FLUID)
C...
C...       NET FLOW RATE OF MOMENTUM WITHIN THE FLUID +
C...
C...       NET PRESSURE FORCE ON THE FLUID
C...
C...       (S*DZ*RHOF*V)
C...                    T
C...
C...       (S*V*RHOF*(V/EPS))     - (S*V*RHOF*(V/EPS))
C...                         AT Z                     AT Z+DZ
C...
C...       - (S*DZ*(DP/DZ)*GC)
C...
C...
C...       DIVISION BY (S*DZ) WITH DZ ----+ 0 GIVES
C...
C...       (RHOF*V)  =
C...               T                                                (5)
C...
C...       - (1/EPS)*(RHOF*(V**2))  - (DP/DZ)*GC
C...                              Z
C...
C...       CHECK ON UNITS
C...
C...       (S*DZ*RHOF*V)
C...                    T
C...
C...       (M**2)*(M)*(KG/M**3)*(M/SEC)*(1/SEC) = KG-M/SEC**2
C...
C...       (S*V*RHOF*(V/EPS))
C...
C...       (M**2)*(M/SEC)*(KG/M**3)*(M/SEC) = KG-M/SEC**2
C...
C...       (S*DZ*(DP/DZ)*GC)
C...
C...       (M**2)*(M)*(N/M**2)*(1/M)*(KG-M/N-SEC**2) = KG-M/SEC**2
C...
C...       (NOTE - 1 N/M**2 = 1 NEWTON/M**2 = 1 PA = 1 PASCAL)
C...
C...  (6)  EQUATION OF STATE FOR THE FLUID
C...       -------------------------------
C...
C...       RHOF = (P*MW)/(R*TF)                                     (6)
C...
C...  (7)  PRESSURE DROP EQUATION
C...       ----------------------
C...
C...       DP/DZ = KL*V + KT*(V**2)                                 (7)
C...
C...       (NOTE - KL AND KT ARE COEFFICIENTS FOR LAMINAR AND
C...               TURBULENT FLOW, RESPECTIVELY)
C...
C...
C...       IF WE ASSUME TURBULENT FLOW, KL = 0 AND
C...
C...       KT = (DP/DZ)/V**2 = (-30000 PA/6 M)*(1/0.3 M/SEC)**2
C...
C...          = -55555.6 (PA-SEC**2)/(M**3)
C...
C...  THE CONSTANTS IN THE EQUATIONS CAN NOW BE EVALUATED FROM THE
C...  GIVEN MODEL PARAMETERS.
C...
C...  EQUATION (1)
C...
C...     H*AV/((1-EPS)*RHOC*CPC) = 60*150/((1-0.95)*6000*400) = 0.075
C...
C...  EQUATION (2)
C...
C...     CPF/(EPS*CVF) = 1.4/0.95 = 1.47
C...
C...     H*AV/(EPS*CVF) = 60*150/(0.95*600/1.4) = 22.105
C...
C...     K*DH/CVF = 12000*70000/(600/1.4) = 19.6E+05
C...
C...  EQUATION (3)
C...
C...     1/EPS = 1/0.95 = 1.053
C...
C...  EQUATION (4)
C...
C...     1/EPS = 1/0.95 = 1.053
C...
C...     K = 400
C...
C...  EQUATION (5)
C...
C...     1/EPS = 1/0.95 = 1.053
C...
C...  EQUATION (6)
C...
C...     R/MW = 8212/40 = 205.3
C...
C...  EQUATION (7)
C...
C...     KL = 0, KT = 55555.6
C...
C...  THE AUXILIARY CONDITIONS FOR EQUATIONS (1) TO (5) ARE
C...
C...     INITIAL CONDITIONS
C...
C...        TC(Z,0) = 573 K
C...
C...        TF(Z,0) = 573 K
C...
C...        RHOF(Z,0) = 5.95 KG/M**3
C...
C...        CA(Z,0) = 0
C...
C...        V(Z,0) = 0
C...
C...     BOUNDARY CONDITIONS
C...
C...        TF(0,T) = 573 K
C...
C...        RHOF(0,T) = 5.95 KG/M**3
C...
C...        V(0,T) = 0.3 M/SEC
C...
C...        CA(0,T) = 20*(1 - EXP(-T/5)) MOLE/M**3
C...
C...  THE FOLLOWING CODING IMPLEMENTS THE PRECEDING MATHEMATICAL MODEL
C...  OVER A 13-POINT NUMERICAL METHOD OF LINES GRID
C...
C...  TEMPORAL COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    TC(13),     R(13),    CA(13),   RTF(13),    RV(13)
     2      /F/ DTCDT(13),  DRDT(13), DCADT(13),DRTFDT(13), DRVDT(13)
C...
C...  SPATIAL COMMON AREA
      COMMON/S/    P1(13),    P2(13),    P3(13),    P4(13),
     1          DP1DZ(13), DP2DZ(13), DP3DZ(13), DP4DZ(13)
C...
C...  VARIABLE COMMON AREA
      COMMON/V/    TF(13),     V(13)
C...
C...  CONSTANT COMMON AREA
      COMMON/C/       TF0,       TC0,        R0,        V0,       CA0,
     1                TFE,                   RE,        VE,       CAE,
     2                  N,        ZL,        DZ
C...
C...  INPUT/OUTPUT DEVICE NUMBERS FOR PRINTING IN DERV
      COMMON/IO/       NI,        NO
C...
C...  GRID PARAMETERS
      N=13
      ZL=6.
      DZ=ZL/FLOAT(N-1)
C...
C...  INITIAL VALUES
      TF0=573.
      TC0=573.
       R0=5.95
       V0=0.3
      CA0=0.
C...
C...  BOUNDARY (ENTERING) VALUES
      TFE=573.
       RE=5.95
       VE=0.3
      CAE=20.
C...
C...  INITIAL CONDITONS
      DO 1 I=1,N
       TF(I)=TF0
       TC(I)=TC0
        R(I)=R0
        V(I)=V0
       CA(I)=CA0
      RTF(I)=R0*TF0
       RV(I)=R0*V0
1     CONTINUE
C...
C...  THE FOLLOWING CALL SO PRINT CAN BE USED TO CHECK THE
C...  INITIAL CONDITIONS SET IN INITAL
C...  CALL PRINT(NI,NO)
C...
C...  INITIALIZE THE DERIVATIVE CALCULATIONS
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
C...
C...  TEMPORAL COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    TC(13),     R(13),    CA(13),   RTF(13),    RV(13)
     2      /F/ DTCDT(13),  DRDT(13), DCADT(13),DRTFDT(13), DRVDT(13)
C...
C...  SPATIAL COMMON AREA
      COMMON/S/    P1(13),    P2(13),    P3(13),    P4(13),
     1          DP1DZ(13), DP2DZ(13), DP3DZ(13), DP4DZ(13)
C...
C...  VARIABLE COMMON AREA
      COMMON/V/    TF(13),     V(13)
C...
C...  CONSTANT COMMON AREA
      COMMON/C/       TF0,       TC0,        R0,        V0,       CA0,
     1                TFE,                   RE,        VE,       CAE,
     2                  N,        ZL,        DZ
C...
C...  INPUT/OUTPUT DEVICE NUMBERS FOR PRINTING IN DERV
      COMMON/IO/       NI,        NO
C...
C...  MODEL CONSTRAINTS
      DO 1 I=1,N
      IF(TF(I).LT.TFE)TF(I)=TFE
      IF(TC(I).LT.TFE)TC(I)=TFE
      IF( R(I).LT. 0.) R(I)=0.
      IF( V(I).LT. 0.) V(I)=0.
      IF(CA(I).LT. 0.)CA(I)=0.
1     CONTINUE
C...
C...  MODEL ALGEBRA
      DO 2 I=1,N
      TF(I)=RTF(I)/R(I)
       V(I)= RV(I)/R(I)
      P1(I)=RTF(I)*V(I)
      P2(I)= RV(I)
      P3(I)= CA(I)*V(I)
      P4(I)= RV(I)*V(I)
2     CONTINUE
C...
C...  MODEL BOUNDARY CONDITIONS
      TF(1)=TFE
       R(1)=RE
       V(1)=VE
      CA(1)=CAE*(1.0-EXP(-T/5.0))
      P1(1)=R(1)*V(1)*TF(1)
      P2(1)=R(1)*V(1)
      P3(1)=V(1)*CA(1)
      P4(1)=R(1)*V(1)**2
C...
C...  MODEL SPATIAL DERIVATIVES
      CALL DSS020(0.,ZL,N,P1,DP1DZ,1.)
      CALL DSS020(0.,ZL,N,P2,DP2DZ,1.)
      CALL DSS020(0.,ZL,N,P3,DP3DZ,1.)
      CALL DSS020(0.,ZL,N,P4,DP4DZ,1.)
C...
C...  MODEL TIME DERIVATIVE AT Z = 0
       DTCDT(1)=0.075*(TF(1)-TC(1))
      DRTFDT(1)=0.
        DRDT(1)=0.
       DCADT(1)=0.
       DRVDT(1)=0.
C...
C...  MODEL TIME DERIVATIVES FOR Z GT 0
      DO 3 I=2,N
       DTCDT(I)=0.075*(TF(I)-TC(I))
             RA=CA(I)*EXP(-7500./TF(I))
      DRTFDT(I)=-1.474*DP1DZ(I)-22.105*(TF(I)-TC(I))+19.6E+05*RA
        DRDT(I)=-1.053*DP2DZ(I)
       DCADT(I)=-1.053*DP3DZ(I)-12000.0*RA
       DRVDT(I)=-1.053*DP4DZ(I)
3     CONTINUE
C...
C...  THE FOLLOWING CALL SO PRINT CAN BE USED TO CHECK THE
C...  CALCULATIONS IN DERV AS THE INTEGRATION PROCEEDS
C...  CALL PRINT(NI,NO)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  TEMPORAL COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    TC(13),     R(13),    CA(13),   RTF(13),    RV(13)
     2      /F/ DTCDT(13),  DRDT(13), DCADT(13),DRTFDT(13), DRVDT(13)
C...
C...  SPATIAL COMMON AREA
      COMMON/S/    P1(13),    P2(13),    P3(13),    P4(13),
     1          DP1DZ(13), DP2DZ(13), DP3DZ(13), DP4DZ(13)
C...
C...  VARIABLE COMMON AREA
      COMMON/V/    TF(13),     V(13)
C...
C...  CONSTANT COMMON AREA
      COMMON/C/       TF0,       TC0,        R0,        V0,       CA0,
     1                TFE,                   RE,        VE,       CAE,
     2                  N,        ZL,        DZ
C...
C...  DIMENSION THE ARRAYS FOR STORING THE NUMERICAL SOLUTION
      DIMENSION TFP(3,101), CAP(3,101), TP(101)
C...
C...  INITIALIZE A COUNTER FOR THE PLOTTED SOLUTION
      DATA IP/0/
C...
C...  MODEL CONSTRAINTS
      DO 4 I=1,N
      IF(TF(I).LT.TFE)TF(I)=TFE
      IF(TC(I).LT.TFE)TC(I)=TFE
      IF( R(I).LT. 0.) R(I)=0.
      IF( V(I).LT. 0.) V(I)=0.
      IF(CA(I).LT. 0.)CA(I)=0.
4     CONTINUE
C...
C...  PRINT THE NUMERICAL SOLUTION IN EVERY 10 CALLS TO PRINT
      IF((IP/10*10-IP).NE.0)GO TO 10
      NS=2
      WRITE(NO,1)T,
     1 (    TC(I),I=1,N,NS),(     R(I),I=1,N,NS),(    CA(I),I=1,N,NS),
     2 (   RTF(I),I=1,N,NS),(    RV(I),I=1,N,NS),( DTCDT(I),I=1,N,NS),
     3 (  DRDT(I),I=1,N,NS),( DCADT(I),I=1,N,NS)
      WRITE(NO,2)
     1 (DRTFDT(I),I=1,N,NS),( DRVDT(I),I=1,N,NS),(    P1(I),I=1,N,NS),
     2 (    P2(I),I=1,N,NS),(    P3(I),I=1,N,NS),(    P4(I),I=1,N,NS),
     3 ( DP1DZ(I),I=1,N,NS),( DP2DZ(I),I=1,N,NS),( DP3DZ(I),I=1,N,NS)
      WRITE(NO,3)
     1 ( DP4DZ(I),I=1,N,NS),(    TF(I),I=1,N,NS),(     V(I),I=1,N,NS)
1     FORMAT(1H1,
     1  5H    T,           F7.2,//,
     2 15H             TC,3E15.3,/,4E15.3,//,
     3 15H              R,3E15.3,/,4E15.3,//,
     4 15H             CA,3E15.3,/,4E15.3,//,
     5 15H            RTF,3E15.3,/,4E15.3,//,
     6 15H             RV,3E15.3,/,4E15.3,//,
     7 15H          DTCDT,3E15.3,/,4E15.3,//,
     8 15H           DRDT,3E15.3,/,4E15.3,//,
     9 15H          DCADT,3E15.3,/,4E15.3,/)
2     FORMAT(
     1 15H         DRTFDT,3E15.3,/,4E15.3,//,
     2 15H          DRVDT,3E15.3,/,4E15.3,//,
     3 15H             P1,3E15.3,/,4E15.3,//,
     4 15H             P2,3E15.3,/,4E15.3,//,
     5 15H             P3,3E15.3,/,4E15.3,//,
     6 15H             P4,3E15.3,/,4E15.3,//,
     7 15H          DP1DZ,3E15.3,/,4E15.3,//,
     8 15H          DP2DZ,3E15.3,/,4E15.3,//,
     9 15H          DP3DZ,3E15.3,/,4E15.3,/)
3     FORMAT(
     1 15H          DP4DZ,3E15.3,/,4E15.3,//,
     2 15H             TF,3E15.3,/,4E15.3,//,
     3 15H              V,3E15.3,/,4E15.3,/)
C...
C...  STORE THE NUMERICAL SOLUTION FOR PLOTTING
10    IP=IP+1
      TFP(1,IP)=TF(5)
      TFP(2,IP)=TF(9)
      TFP(3,IP)=TF(13)
      CAP(1,IP)=CA(5)
      CAP(2,IP)=CA(9)
      CAP(3,IP)=CA(13)
      TP(IP)=T
C...
C...  PLOT THE NUMERICAL SOLUTION
      IF(IP.LT.101)RETURN
      CALL SPLOTS(3,IP,TP,TFP)
      WRITE(NO,11)
11    FORMAT(1H ,/,
     1 12H 1 - TF(2,T),/,
     2 12H 2 - TF(4,T),/,
     3 12H 3 - TF(6,T))
      CALL SPLOTS(3,IP,TP,CAP)
      WRITE(NO,12)
12    FORMAT(1H ,/,
     1 12H 1 - CA(2,T),/,
     2 12H 2 - CA(4,T),/,
     3 12H 3 - CA(6,T))
      RETURN
      END

