*DECK PRO2P5
      SUBROUTINE INITAL
C...
C...  SUBROUTINE INITAL IS CALLED BY MAIN PROGRAM PRO1P3 TO DEFINE THE
C...  INITIAL CONDITION IN THE NUMOL SOLUTION OF EQUATIONS (2.17) TO
C...  (2.20)
C...
      COMMON/T/   TIME,    NSTOP,    NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /R/  R(51),   TR(51),  TRR(51),
     4              R0,       T0,       TA,      A,      B,     RG,
     5               E,       DR
     6      /I/      N
      REAL           K,       K0
C...
C...  SET THE PROBLEM CONSTANTS (ENERGY IS IN CALORIES RATHER THAN
C...  JOULES)
      T0=25.
      R0=1.0
      E=30800.0
      RG=1.987
      RHO=1.71
      CP=0.300
      K=2.4E-04
      DH=1276.
      K0=1.26E+11
C...
C...  SET THE SURFACE TEMPERATURE FOR BOUNDARY CONDITION (2.20)
      IF(NORUN.EQ.1)TA=125.
      IF(NORUN.EQ.2)TA=158.
C...
C...  COMPUTE THE THERMAL DIFFUSIVITY, COEFFICIENT OF REACTION RATE TERM
      A=K/(RHO*CP)
      B=(DH*K0/(RHO*CP))
C...
C...  GRID SPACING
      N=51
      DR=R0/FLOAT(N-1)
C...
C...  INITIAL CONDITION (2.18)
      DO 1 I=1,N
         R(I)=FLOAT(I-1)*DR
         T(I)=T0
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV IS CALLED BY INTEGRATOR RKF45 TO DEFINE THE
C...  ODES IN THE NUMOL SOLUTION OF EQUATIONS (2.17) TO (2.20)
C...
      COMMON/T/   TIME,    NSTOP,    NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /R/  R(51),   TR(51),  TRR(51),
     4              R0,       T0,       TA,      A,      B,     RG,
     5               E,       DR
     6      /I/      N
      REAL           K,       K0
C...
C...  BOUNDARY CONDITION (2.20)
       T(N)=TA
      TT(N)=0.
C...
C...  DERIVATIVE TR
      CALL DSS002(0.,R0,N,T,TR)
C...
C...  BOUNDARY CONDITION (2.19)
      TR(1)=0.
C...
C...  DERIVATIVE TRR
      CALL DSS002(0.,R0,N,TR,TRR)
C...
C...  EQUATION (2.17)
      DO 1 I=1,N-1
C...
C...     R NE 0
         IF(I.NE.1)THEN
            TT(I)=A*(TRR(I)+(2./R(I))*TR(I))
     1           +B*EXP(-E/(RG*(T(I)+273.16)))

C...     R = 0
         ELSE IF(I.EQ.1)THEN
            TT(I)=3.*A*TRR(I)
     1              +B*EXP(-E/(RG*(T(I)+273.16)))
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  SUBROUTINE PRINT IS CALLED BY MAIN PROGRAM PRO1P3 TO PRINT THE
C...  NUMOL SOLUTION OF EQUATIONS (2.17) TO (2.20)
C...
      COMMON/T/   TIME,    NSTOP,    NORUN
     1      /Y/  T(51)
     2      /F/ TT(51)
     3      /R/  R(51),   TR(51),  TRR(51),
     4              R0,       T0,       TA,      A,      B,     RG,
     5               E,       DR
     6      /I/      N
      REAL           K,       K0
C...
C...  PRINT THE NUMOL SOLUTION TO EQUATIONS (2.17) TO (2.20)
      WRITE(NO,2)TIME,(T(I),I=1,N,10)
2     FORMAT(' TIME = ',F6.1,/,15X,'  R=0',5X,'R=0.2',5X,'R=0.4',5X,
     1                             'R=0.6',5X,'R=0.8',5X,'  R=1',/,
     2                         '    T(R,T)',6F10.3,/)
      RETURN
      END
