*DECK SEP34
      SUBROUTINE INITAL
C...
C...  RELATIVE EFFECTS OF SURFACE MASS TRANSFER RESISTANCE, INTERNAL
C...  DIFFUSION AND INTERNAL REACTION ON THE PERFORMANCE OF A SPHERICAL
C...  CATALYST PARTICLE
C...
C...  CONSIDER THE FOLLOWING MODEL FOR DIFFUSION AND REACTION IN A
C...  SPHERICAL CATALYST PARTICLE CONTAINING AN IMMOBILIZED ENZYME:
C...
C...  C  = D*(C   + (2/R)*C ) - K*C
C...   T       RR          R
C...
C...  C(R,0) = C0, C (0,T) = 0, -D*C (R0,T) = KY*(C(R0,T) - CS)
C...                R               R
C...
C...  WHERE
C...
C...  C        CONCENTRATION OF THE SUBSTRATE (REACTANT)
C...
C...  R        RADIAL POSITION WITHIN THE PARTICLE
C...
C...  T        TIME
C...
C...  D        DIFFUSIVITY WITHIN THE PARTICLE
C...
C...  K        FIRST-ORDER REACTION RATE CONSTANT
C...
C...  R0       RADIUS OF THE PARTICLE
C...
C...  KY       MASS TRANSFER COEFFICIENT
C...
C...  CS       SUBSTRATE CONCENTRATION IN THE LIQUID SURROUNDING THE
C...           PARTICLE
C...
C...  C0       INITIAL SUBSTRATE CONCENTRATION
C...
C...  (A)  EXPRESS THE PRECEDING SET OF EQUATIONS IN TERMS OF THE
C...       DIMENSIONLESS VARIABLES
C...
C...       CD = C/CS, TD = T*D/R0**2, RD = R/R0
C...
C...  (B)  COMPUTE A SOLUTION TO THE DIMENSIONLESS EQUATIONS FOR
C...       C0 = 0, CS = 0.001, K = D = 1, R0 = 0.1, KY = 0.1.  PLOT
C...       CD(RD,TD) VS RD FOR A SERIES OF DIMENSIONLESS TIMES, TD.
C...       USE SUBROUTINE DSS004 FOR FIRST-ORDER DERIVATIVES AND DSS042
C...       FOR SECOND-ORDER DERIVATIVES IN RD.
C...
C...  THE SYSTEM WHEN WRITTEN IN TERMS OF THE DIMENSIONLESS VARIABLES
C...  IS
C...
C...  CD  = (CD     + (2/RD)*CD  ) - K*(R0**2)/D*CD
C...    TD     RDRD            RD
C...
C...  CD(RD,0) = C0/CS, CD  (0,TD) = 0
C...                      RD
C...
C...  -CD  (1,T) = KY*R0/D*(CD(1,T) - 1)
C...     RD
C...
C...  NOTE THE IMPORTANT DIMENSIONLESS GROUPS K*(R0**2)/D AND KY*R0/D,
C...  WHICH WILL BE DESIGNATED A1 AND A2, RESPECTIVELY, IN THE FOLLOW-
C...  ING CODING.  ALSO, THE DIMENSIONLESS DESIGNATION, D, FOR THE
C...  VARIABLES WILL BE DROPPED TO SIMPLIFY THE CODING, E.G., CD WILL
C...  BE REPLACED WITH C, WITH THE UNDERSTANDING THAT THIS IS A DI-
C...  MENSIONLESS CONCENTRATION.  THE SAME APPLIES TO TD (NOW T) AND
C...  RD (NOW R).  AN ADVANTAGE OF USING DIMENSIONLESS VARIABLES IS
C...  THE RANGE OVER WHICH THEY VARY, E.G., 0 LE R LE 1 REGARDLESS OF
C...  THE VALUE OF R0.
C...
C...  FINALLY, NOTE ONE OTHER IMPORTANT CHARACTERISTIC OF THIS SYSTEM:
C...  THE CONCENTRATION CS USED TO DEFINE THE DIMENSIONLESS CONCENTRA-
C...  TION, DOES NOT APPEAR IN THE FINAL DIMENSIONLESS EQUATIONS.  THIS
C...  IS BECAUSE THE PROBLEM IS ENTIRELY LINEAR.  IF THE PROBLEM HAD
C...  BEEN NONLINEAR, CS WOULD APPEAR SOMEWHERE IN THE EQUATIONS, I.E.,
C...  A FUNDAMENTAL PROPERTY OF NONLINEAR SYSTEMS IS THE DEPENDENCE ON
C...  THE ABSOLUTE VALUES OF THE STATE VARIABLES.  TO GIVE A SPECIFIC
C...  EXAMPLE, CONSIDER THE FIRST EQUATION IF THE REACTION HAD BEEN
C...  SECOND ORDER, SO THAT THE REACTION TERM WAS -K*(C**2).  NOW IF
C...  THE SAME DIMENSIONLESS VARIABLES ARE USED, CS WILL APPEAR IN THE
C...  DIMENSIONLESS PDE.
C...
C...  DSS/2 COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(11)
     2      /F/    CT(11)
     3      /S/    CR(11),   CRR(11)
     4      /C/        A1,        A2,       N11,        IP
      REAL              K,        KY
C...
C...  SET THE MODEL PARAMETERS
      C0=0.
      CS=0.001
      R0=0.1
C...
C...  SLOW MASS TRANSFER, SLOW REACTION, FAST DIFFUSION
      IF(NORUN.EQ.1)THEN
      KY=0.1
      K =1.0
      D =1.0
C...
C...  FAST MASS TRANSFER, FAST REACTION, SLOW DIFFUSION
      ELSE IF(NORUN.EQ.2)THEN
      KY=5.0
      K =50.0
      D =0.1
      END IF
C...
C...  DIMENSIONLESS GROUPS
      A1=K*(R0**2)/D
      A2=KY*R0/D
C...
C...  INITIAL CONDITION
      N11=11
      DO 1 I=1,N11
      C(I)=0.
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(11)
     2      /F/    CT(11)
     3      /S/    CR(11),   CRR(11)
     4      /C/        A1,        A2,       N11,        IP
      REAL              K,        KY
C...
C...  FIRST-ORDER SPATIAL DERIVATIVE
      CALL DSS004(0.,1.,N11,C,CR)
C...
C...  BOUNDARY CONDITIONS
      CR(1)=0.
      CR(N11)=-A2*(C(N11)-1.0)
C...
C...  SECOND-ORDER SPATIAL DERIVATIVE
      NL=2
      NU=2
      CALL DSS042(0.,1.,N11,C,CR,CRR,NL,NU)
C...
C...  PDE, R NE 0
      DO 1 I=2,N11
      R=FLOAT(I-1)/10.
      CT(I)=(CRR(I)+(2.0/R)*CR(I))-A1*C(I)
1     CONTINUE
C...
C...  PDE, R = 0
      CT(1)=3.0*CRR(1)-A1*C(1)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(11)
     2      /F/    CT(11)
     3      /S/    CR(11),   CRR(11)
     4      /C/        A1,        A2,       N11,        IP
      REAL              K,        KY
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION RP(11), CP(2,11)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,1)
1        FORMAT(9X,'T','   (0.0,T)','   (0.2,T)','   (0.4,T)',
     1                 '   (0.6,T)','   (0.8,T)','   (1.0,T)')
      END IF
C...
C...  PRINT THE SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
         WRITE(NO,2)T,(  C(I),I=1,N11,2),
     1                ( CR(I),I=1,N11,2),
     2                (CRR(I),I=1,N11,2),
     3                ( CT(I),I=1,N11,2)
2     FORMAT(F10.2,/,'    C(R,T)',6F10.4,/,
     1               '   CR(R,T)',6F10.4,/,
     2               '  CRR(R,T)',6F10.4,/,
     3               '   CT(R,T)',6F10.4,/)
      END IF
C...
C...  STORE THE CONCENTRATION PROFILE AT THE END OF THE RUN FOR
C...  PLOTTING
      IF(IP.EQ.101)THEN
         DO 3 I=1,N11
         CP(NORUN,I)=C(I)
         RP(I)=FLOAT(I-1)/10.0
3        CONTINUE
C...
C...     STORE THE DAMKOHLER AND BIOT NUMBERS FOR SUBSEQUENT PRINTING
         IF(NORUN.EQ.1)D1=A1
         IF(NORUN.EQ.1)B1=A2
         IF(NORUN.EQ.2)D2=A1
         IF(NORUN.EQ.2)B2=A2
      END IF
C...
C...  PLOT THE CONCENTRATION PROFILES AT THE END OF THE SECOND RUN
      IF((IP.EQ.101).AND.(NORUN.EQ.2))THEN
      CALL SPLOTS(NORUN,N11,RP,CP)
      WRITE(NO,11)D1,D2,B1,B2
11    FORMAT(' C(R,T) VS R',
     1       '   1 - MASS TRANSFER AND REACTION LIMITED',
     2       '   2 - DIFFUSION LIMITED',/,
     3       ' DAMKOHLER NUMBERS = ',2E10.3,
     4       '   BIOT NUMBERS = '   ,2E10.3)
      END IF
      RETURN
      END
