*DECK REACT28
      SUBROUTINE INITAL
C...
C...  METHOD OF LINES ANALYSIS OF THE OUTFLOW BOUNDARY CONDITION FOR
C...  A ONE-DIMENSIONAL SYSTEM WITH CONVECTION AND AXIAL DISPERSION
C...
C...  THE FOLLOWING MATHEMATICAL MODEL HAS BEEN PROPOSED FOR A TUBULAR
C...  REACTOR WITH AXIAL DISPERSION
C...
C...     C  = -V*C  + D*C   - K*C                               (1)
C...      T       Z      ZZ
C...
C...  WHERE
C...
C...     C       CONCENTRATION OF REACTANT
C...
C...     T       TIME
C...
C...     Z       AXIAL POSITION ALONG THE REACTOR
C...
C...     V       AVERAGE VELOCITY THROUGH THE REACTOR
C...
C...     D       AXIAL DISPERSION COEFFICIENT
C...
C...     K       REACTION VELOCITY CONSTANT
C...
C...  EQUATION (1) CAN BE APPLIED TO OTHER PHYSICAL SYSTEMS.  FOR EXAMPLE,
C...  IF WE CONSIDER A HEAT TRANSFER SYSTEM WITH FLUID TEMPERATURE TM,
C...  THEN AN ENERGY BALANCE GIVES
C...
C...     TM  = -V*TM  + D*TM   - UA*TM
C...       T        Z       ZZ
C...
C...  WHERE UA IS THE PRODUCT OF THE HEAT TRANSFER AREA AND HEAT TRANSFER
C...  COEFFICIENT.  THE TEMPERATURE DIFFERENCE IN THIS CASE IS (TM - 0)
C...  WHICH IS REASONABLE FOR A LINEAR SYSTEM SINCE THE ENERGY BALANCE
C...  CAN BE WRITTEN RELATIVE TO A ZERO TEMPERATURE (I.E., THE ABSOLUTE
C...  MAGNITUDE OF THE DATUM TEMPERATURE IS IMMATERIAL FOR A LINEAR SYS-
C...  TEM).
C...
C...  THE INITIAL AND BOUNDARY CONDITIONS FOR EQUATION (1) ARE
C...
C...     C(Z,0) = C0, C(0,T) = CIN                           (2)(3)
C...
C...     CASE 1: DIRICHLET BOUNDARY CONDITION AT Z = ZL
C...
C...             C(ZL,T) = COUT                                 (4)
C...
C...     CASE 2: NEUMANN BOUNDARY CONDITION AT Z = ZL
C...
C...             C (ZL,T) = 0                                   (5)
C...              Z
C...
C...     CASE 3: BOUNDARY CONDITION OF THE THIRD TYPE AT Z = ZL
C...
C...             C (ZL,T) = -V*C (ZL,T) - K*C(ZL,T)             (6)
C...              T             Z
C...     
C...  IF DIMENSIONLESS VARIABLES ARE DEFINED AS
C...
C...     U=C/CIN, Z = Z/ZL, T = T*D/ZL**2
C...
C...  EQUATIONS (1) CAN BE WRITTEN AS
C...
C...     CIN*(D/ZL**2)*U  = - V*(CIN*/ZL)*U  + D*(CIN/ZL**2)*U   
C...                    T                  Z                  ZZ  
C...                        - K*CIN*U
C...  OR
C...
C...     U   = -PE*U  + U   - DA*U                              (7)
C...      T         Z    ZZ
C...
C...  WHERE
C...
C...     PE       V*ZL/D (PECLET NUMBER)
C...
C...     DA       K*(ZL**2)/D (DAMKOHLER NUMBER)
C...
C...  SIMILARLY, EQUATIONS (2), (3), (4), (5) AND (6) BECOME
C...
C...     U(Z,0) = C0/CIN = U0, U(0,T) = 1                    (8)(9)
C...
C...     CASE 1:
C...
C...     U(1,T) = COUT/CIN = UL                                (10)
C...
C...     CASE 2:
C...
C...     U (1,T) = 0                                           (11)
C...      Z
C...
C...     CASE 3:
C...
C...     U (1,T) = -PE*U (1,T) - DA*U(1,T)                     (12)
C...      T             Z
C...
C...  WE NOW CONSIDER THE SOLUTION OF EQUATIONS (7) TO (12) FOR 
C...  A SERIES OF PARAMETRIC VALUES OF PE, DA AND U0.  ALSO, THE
C...  RESULTS AT STEADY STATE ARE COMPARED WITH THE ANALYTICAL
C...  SOLUTIONS DERVIVED BELOW.  EQUATION (7) AT STEADY STATE IS 
C...
C...     0 = -PE*U  + U   - DA*U                               (13)
C...              Z    ZZ
C...  
C...  IF A STEADY STATE SOLUTION TO EQUATION (13) OF THE FORM
C...
C...     U(Z) = C*EXP(M*Z)                                     (14)
C...
C...  IS ASSUMED, SUBSTITUTION OF EQUATION (14) IN EQUATION (13) 
C...  GIVES
C...
C...     0 = -PE*C*M*EXP(M*Z) + C*(M**2)*EXP(M*Z) - DA*C*EXP(M*Z)
C...
C...  OR
C...
C...     M**2 - PE*M - DA = 0                                  (15)
C...
C...
C...  EQUATION (15) IS THE CHARACTERISTIC EQUATION FOR EQUATION (13),
C...  WITH EIGENVALUES 
C...
C...     M1,M2 = (PE + (OR -)(PE**2 - 4*(-DA))**0.5)/2         (16)
C...
C...  THEN, THE SOLUTION TO EQUATION (13) (IF PE NE 0 SO M1 AND M2
C...  ARE DISTINCT) IS
C...
C...     U(Z) = C1*EXP(M1*Z) + C2*EXP(M2*Z)                    (17)
C...
C...  WE NOW CONSIDER APPLICATION OF THE BOUNDARY CONDITIONS FOR THE
C...  THREE CASES.
C...
C...  CASE 1: (BOUNDARY CONDITIONS (9) AND (10))
C...
C...     U(0) = 1 = C1 + C2
C...
C...     U(1) = UL = C1*EXP(M1) + C2*EXP(M2)
C...
C...  WHICH CAN BE SOLVED FOR C1 AND BC2
C...
C...     C1 =  (UL - EXP(M2))/(EXP(M1) - EXP(M2))              (18)
C...
C...     C2 = -(UL - EXP(M1))/(EXP(M1) - EXP(M2))              (19)
C...
C...  CASE 2: (BOUNDARY CONDITIONS (9) AND (11) - SAME AS CASE 3, 
C...           AT STEADY STATE, WITH DA = 0 IN EQUATION (12))
C...
C...  CASE 3: (BOUNDARY CONDITIONS (9) AND (12))
C...
C...     U(0) = 1 = C1 + C2
C...
C...     0 = -PE*U (1) - DA*U(1) = 
C...              Z
C...
C...         - PE*(C1*M1*EXP(M1) + C2*M2*EXP(M2))
C...
C...         - DA*(C1*EXP(M1) + C2*EXP(M2))
C...
C...       = -(PE*M1 + DA)*EXP(M1)*C1 
C...
C...         -(PE*M2 + DA)*EXP(M2)*C2
C...
C...       = A*C1 + B*C2
C...
C...  WHICH CAN BE SOLVED FOR C1 AND C2
C...
C...     C1 = -B/(A - B)                                       (20)
C...
C...     C2 =  A/(A - B)                                       (21)
C...
C...     A = -(PE*M1 + DA)*EXP(M1)  
C...
C...     B =  (PE*M2 + DA)*EXP(M2)
C...
C...  THESE THREE CASES ARE COMPUTED IN SUBROUTINES UA1, UA2 AND UA3
C...  FOR COMPARISON WITH THE NUMERICAL SOLUTION AT STEADY STATE.  
C...
C...  FINALLY, FOR DA = 0, EQUATION (16) GIVES M1 = PE, M2 = 0.  
C...  THEN, FOR CASE 1, EQUATIONS (17), (18) AND (19) REDUCE TO
C...
C...     U(Z) = C1*EXP(PE*Z) + C2                              (22)
C...
C...     C1 =  (UL - 1)/(EXP(PE) - 1)                          (23)
C...
C...     C2 = -(UL - EXP(PE))/(EXP(PE) - 1)                    (24)
C...
C...  COMBINATION OF EQUATIONS (22), (23) AND (24) GIVES
C...
C...     U(Z) = ((UL - 1)*EXP(PE*Z) - (UL - EXP(PE)))/(EXP(PE) - 1)
C...
C...  OR
C...
C...     (U(Z) - 1)/(UL - 1) = (EXP(PE*Z) - 1)/(EXP(PE) - 1))  (25)
C...
C...  WHICH IS THE SOLUTION REPORTED BY PATANKAR, S. V., NUMERICAL
C...  HEAT TRANSFER AND FLUID FLOW, HEMISPHERE PUBLISHING CORPORATION,
C...  NY, 1980, EQ (5.17), P85.  NOTE THAT EQUATION (25) IS INDETER-
C...  MINATE FOR PE = 0.  IF L'HOSPITAL'S RULE IS APPLIED TO THE RHS
C...  OF EQUATION (25)
C...
C...     LIM (PE --> 0) = Z*EXP(Z*PE)/EXP(PE) = Z
C...
C...  AND EQUATION (25) BECOMES
C...
C...     (U(Z) - 1)/(UL - 1) = Z                               (26)
C...
C...  THIS RESULT IS EXPECTED SINCE EQUATION (7) FOR PE = DA = 0 REDUCES 
C...  TO THE ONE-DIMENSIONAL HEAT CONDUCTION EQUATION (FOURIER'S SECOND 
C...  LAW) (ALSO THE ONE-DIMENSIONAL DIFFUSION EQUATION (FICK'S SECOND 
C...  LAW)), FOR WHICH THE STEADY STATE SOLUTION IS LINEAR FROM U(0) TO 
C...  U(1)).
C...
C...  FOR CASES 2 AND 3 WITH DA = 0, EQUATIONS (17), (20) AND (21) REDUCE 
C...  TO (AGAIN WITH M1 = PE, M2 = 0 FROM EQUATION (16))
C...
C...     U(Z) = 1                                              (27)
C...
C...     C1 = 0, C2 = 1
C...
C...  WHICH IS NOT A PARTICULARLY INTERESTING SOLUTION.
C...
C...  CASES 2 AND 3 PROVIDE MORE INTERESTING STEADY STATE SOLUTIONS FOR
C...  DA NE 0 AS DEMONSTRATED BY THE FOLLOWING PROGRAMMING.
C...
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  CASE 1: (EQUATIONS (7), (8), (9) AND (10))
      IF((NORUN.GE.1).AND.(NORUN.LE.3))CALL INIT1
C...
C...  CASE 2: (EQUATIONS (7), (8), (9) AND (11))
      IF((NORUN.GE.4).AND.(NORUN.LE.6))CALL INIT2
C...
C...  CASE 3: (EQUATIONS (7), (8), (9) AND (12))
      IF((NORUN.GE.7).AND.(NORUN.LE.9))CALL INIT3
      RETURN
      END
      SUBROUTINE INIT1
C...
C...  COMMON AREA
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  SET THE MODEL PARAMETERS
      IF(NORUN.EQ.1)THEN
         DA=0.
         PE=0.
         UL=2.0
      ELSE IF(NORUN.EQ.2)THEN
         DA=0.
         PE=10.0
         UL=2.0
      ELSE IF(NORUN.EQ.3)THEN
         DA=0.
         PE=25.0
         UL=2.0
      END IF
C...
C...  MODEL INITIAL CONDITION
      DO 1 I=1,NZ
      U(I)=1.0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT2
C...
C...  COMMON AREA
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  SET THE MODEL PARAMETERS
      IF(NORUN.EQ.4)THEN
         DA=10.0
         PE=1.0
      ELSE IF(NORUN.EQ.5)THEN
         DA=10.0
         PE=10.0
      ELSE IF(NORUN.EQ.6)THEN
         DA=10.0
         PE=25.0
      END IF
C...
C...  MODEL INITIAL CONDITION
      DO 1 I=1,NZ
      U(I)=1.0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE INIT3
C...
C...  COMMON AREA
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  SET THE MODEL PARAMETERS
      IF(NORUN.EQ.7)THEN
         DA=10.0
         PE=1.0
      ELSE IF(NORUN.EQ.8)THEN
         DA=10.0
         PE=10.0
      ELSE IF(NORUN.EQ.9)THEN
         DA=10.0
         PE=25.0
      END IF
C...
C...  MODEL INITIAL CONDITION
      DO 1 I=1,NZ
      U(I)=1.0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  CASE 1: (EQUATIONS (7), (8), (9) AND (10))
      IF((NORUN.GE.1).AND.(NORUN.LE.3))CALL DERV1
C...
C...  CASE 2: (EQUATIONS (7), (8), (9) AND (11))
      IF((NORUN.GE.4).AND.(NORUN.LE.6))CALL DERV2
C...
C...  CASE 3: (EQUATIONS (7), (8), (9) AND (12))
      IF((NORUN.GE.7).AND.(NORUN.LE.9))CALL DERV3
      RETURN
      END
      SUBROUTINE DERV1
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=1
      U(1)=1.0
C...
C...  BOUNDARY CONDITION AT Z = 1
      NU=1
      U(NZ)=UL
C...
C...  U
C...   Z
      CALL DSS020(0.,1.0,NZ,U,UZ,1.0)
C...
C...  U
C...   ZZ
      CALL DSS044(0.,1.0,NZ,U,UZ,UZZ,NL,NU)
C...
C...  PDE
      DO 1 I=1,NZ
C...
C...     Z = 0
         IF(I.EQ.1)THEN
            UT(I)=0.
C...
C...     Z = 1
         ELSE IF(I.EQ.NZ)THEN
            UT(I)=0.
C...
C...     INTERIOR Z
         ELSE
            UT(I)=-PE*UZ(I)+UZZ(I)-DA*U(I)         
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=1
      U(1)=1.0
C...
C...  U
C...   Z
      CALL DSS020(0.,1.0,NZ,U,UZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 1 
      NU=2
      UZ(NZ)=0.
C...
C...  U
C...   ZZ
      CALL DSS044(0.,1.0,NZ,U,UZ,UZZ,NL,NU)
C...
C...  PDE
      DO 1 I=1,NZ
C...
C...     Z = 0
         IF(I.EQ.1)THEN
            UT(I)=0.
C...
C...     Z NE 0
         ELSE
            UT(I)=-PE*UZ(I)+UZZ(I)-DA*U(I)         
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  BOUNDARY CONDITION AT Z = 0
      NL=1
      U(1)=1.0
C...
C...  U
C...   Z
      CALL DSS020(0.,1.0,NZ,U,UZ,1.0)
C...
C...  BOUNDARY CONDITION AT Z = 1.  UZ(NZ) FROM THE PREVIOUS CALL TO 
C...  DSS020 IS USED AS THE BOUNDARY CONDITION FOR THE FOLLOWING CALL
C...  TO DSS044.  AN ALTERNATIVE APPROACH IS TO SPECIFY A DIRICHLET
C...  BOUNDARY CONDITION (NU=1) AND USE U(NZ) FROM COMMON/Y/ AS A
C...  BOUNDARY VALUE IN THE FOLLOWING CALL TO DSS044
      NU=2
C...
C...  U
C...   ZZ
      CALL DSS044(0.,1.0,NZ,U,UZ,UZZ,NL,NU)
C...
C...  PDE
      DO 1 I=1,NZ
C...
C...     Z = 0
         IF(I.EQ.1)THEN
            UT(I)=0.
C...
C...     Z = 1
         ELSE IF(I.EQ.NZ)THEN
            UT(I)=-PE*UZ(I)-DA*U(I)
C...
C...     INTERIOR Z
         ELSE
            UT(I)=-PE*UZ(I)+UZZ(I)-DA*U(I)         
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  CASE 1: (EQUATIONS (7), (8), (9) AND (10))
      IF((NORUN.GE.1).AND.(NORUN.LE.3))CALL PRINT1(NI,NO)
C...
C...  CASE 2: (EQUATIONS (7), (8), (9) AND (11))
      IF((NORUN.GE.4).AND.(NORUN.LE.6))CALL PRINT2(NI,NO)
C...
C...  CASE 3: (EQUATIONS (7), (8), (9) AND (12))
      IF((NORUN.GE.7).AND.(NORUN.LE.9))CALL PRINT3(NI,NO)
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION ZP(NZ), UP(2,NZ)
C...
C...  MONITOR THE PROGRESS OF THE SOLUTION
      IP=IP+1
      WRITE(*,*)' NORUN = ',NORUN,'  IP = ',IP
C...
C...  PRINT THE NUMERICAL SOLUTION
      NZI=1
      WRITE(NO,1)T,(U(I),I=1,NZ,NZI)
1     FORMAT(/,' T = ',F5.1,/,'    U(Z,T)',/,F10.4,/,(5F10.4))
C...
C...  STORE THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS FOR 
C...  PLOTTING
      IF(IP.EQ.6)THEN
      DO 2 I=1,NZ
         Z=FLOAT(I-1)/FLOAT(NZ-1)
         ZP(I)=Z
         UP(1,I)=U(I)
         UP(2,I)=UA1(PE,DA,UL,Z)
2     CONTINUE
      END IF
C...
C...  PRINT THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS AND THE
C...  DIFFERENCE
      IF(IP.LT.6)RETURN
      WRITE(NO,3)(I,ZP(I),UP(1,I),UP(2,I),UP(1,I)-UP(2,I),I=1,NZ)
3     FORMAT(//,4X,'I',9X,'Z',2X,' U (NUM)',2X,' U (ANAL)',
     +      2X,' U (DIFF)',/,(I5,4F10.4))
C...
C...  PLOT THE SPATIAL PROFILES, U(Z,INF) VS Z, AT STEADY STATE
      CALL TPLOTS(2,NZ,ZP,UP)
      WRITE(NO,11)
11    FORMAT(1H ,/,'   U(Z,T) VS Z, 1 - NUMERICAL, 2 - ANALYTICAL')
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION ZP(NZ), UP(2,NZ)
C...
C...  MONITOR THE PROGRESS OF THE SOLUTION
      IP=IP+1
      WRITE(*,*)' NORUN = ',NORUN,'  IP = ',IP
C...
C...  PRINT THE NUMERICAL SOLUTION
      NZI=1
      WRITE(NO,1)T,(U(I),I=1,NZ,NZI)
1     FORMAT(/,' T = ',F5.1,/,'    U(Z,T)',/,F10.4,/,(5F10.4))
C...
C...  STORE THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS FOR 
C...  PLOTTING
      IF(IP.EQ.6)THEN
      DO 2 I=1,NZ
         Z=FLOAT(I-1)/FLOAT(NZ-1)
         ZP(I)=Z
         UP(1,I)=U(I)
         UP(2,I)=UA2(PE,DA,UL,Z)
2     CONTINUE
      END IF
C...
C...  PRINT THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS AND THE
C...  DIFFERENCE
      IF(IP.LT.6)RETURN
      WRITE(NO,3)(I,ZP(I),UP(1,I),UP(2,I),UP(1,I)-UP(2,I),I=1,NZ)
3     FORMAT(//,4X,'I',9X,'Z',2X,' U (NUM)',2X,' U (ANAL)',
     +      2X,' U (DIFF)',/,(I5,4F10.4))
C...
C...  PLOT THE SPATIAL PROFILES, U(Z,INF) VS Z, AT STEADY STATE
      CALL TPLOTS(2,NZ,ZP,UP)
      WRITE(NO,11)
11    FORMAT(1H ,/,'   U(Z,T) VS Z, 1 - NUMERICAL, 2 - ANALYTICAL')
      RETURN
      END
      SUBROUTINE PRINT3(NI,NO)
      PARAMETER (NZ=41)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(NZ)
     2      /F/    UT(NZ)
     3      /S/    UZ(NZ),   UZZ(NZ)
     4      /C/        PE,        DA,        UL
     5      /I/        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION ZP(NZ), UP(2,NZ)
C...
C...  MONITOR THE PROGRESS OF THE SOLUTION
      IP=IP+1
      WRITE(*,*)' NORUN = ',NORUN,'  IP = ',IP
C...
C...  PRINT THE NUMERICAL SOLUTION
      NZI=1
      WRITE(NO,1)T,(U(I),I=1,NZ,NZI)
1     FORMAT(/,' T = ',F5.1,/,'    U(Z,T)',/,F10.4,/,(5F10.4))
C...
C...  STORE THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS FOR 
C...  PLOTTING
      IF(IP.EQ.6)THEN
      DO 2 I=1,NZ
         Z=FLOAT(I-1)/FLOAT(NZ-1)
         ZP(I)=Z
         UP(1,I)=U(I)
         UP(2,I)=UA3(PE,DA,UL,Z)
2     CONTINUE
      END IF
C...
C...  PRINT THE STEADY STATE NUMERICAL AND ANALYTICAL SOLUTIONS AND THE
C...  DIFFERENCE
      IF(IP.LT.6)RETURN
      WRITE(NO,3)(I,ZP(I),UP(1,I),UP(2,I),UP(1,I)-UP(2,I),I=1,NZ)
3     FORMAT(//,4X,'I',9X,'Z',2X,' U (NUM)',2X,' U (ANAL)',
     +      2X,' U (DIFF)',/,(I5,4F10.4))
C...
C...  PLOT THE SPATIAL PROFILES, U(Z,INF) VS Z, AT STEADY STATE
      CALL TPLOTS(2,NZ,ZP,UP)
      WRITE(NO,11)
11    FORMAT(1H ,/,'   U(Z,T) VS Z, 1 - NUMERICAL, 2 - ANALYTICAL')
      RETURN
      END
      FUNCTION UA1(PE,DA,UL,Z)
C...
C...  FUNCTION UA1 COMPUTES THE STEADY STATE SOLUTION FOR CASE 1
C...
      REAL M1, M2
C...
C...  EQUATION (26)
      IF((PE.EQ.0.).AND.(DA.EQ.0.))THEN
         UA1=1.0+(UL-1.0)*Z
         RETURN
C...  
C...  EQUATION (16)
      ELSE
         M1=(PE+SQRT(PE**2+4.0*DA))/2.0
         M2=(PE-SQRT(PE**2+4.0*DA))/2.0
C...
C...  EQUATIONS (18) AND (19) (NOTE - FOR LARGE PE, E.G., PE = 100,
C...  THESE EXPONENTIALS OVERFLOW, BUT GENERALLY, THE NUMERICAL SOLU-
C...  TION CAN BE COMPUTED)
         EXPM1=EXP(M1)
         EXPM2=EXP(M2)
         C1= (UL-EXPM2)/(EXPM1-EXPM2)
         C2=-(UL-EXPM1)/(EXPM1-EXPM2)
C...
C...  EQUATION (17)
         UA1=C1*EXP(M1*Z)+C2*EXP(M2*Z)
         RETURN
      END IF
      END
      FUNCTION UA2(PE,DA,UL,Z)
C...
C...  FUNCTION UA2 COMPUTES THE STEADY STATE SOLUTION FOR CASE 2
C...
      REAL M1, M2
C...
C...  EQUATION (27)
      IF(DA.EQ.0.)THEN
         UA2=1.0
         RETURN
C...  
C...  EQUATION (16)
      ELSE
         M1=(PE+SQRT(PE**2+4.0*DA))/2.0
         M2=(PE-SQRT(PE**2+4.0*DA))/2.0
C...
C...  EQUATIONS (20) AND (21) FOR DA = 0 AS REQUIRED BY BOUNDARY CONDI-
C...  TION (11) RATHER THAN (12) (NOTE - FOR LARGE PE, E.G., PE = 100,
C...  THESE EXPONENTIALS OVERFLOW, BUT GENERALLY, THE NUMERICAL SOLU-
C...  TION CAN BE COMPUTED)
         EXPM1=EXP(M1)
         EXPM2=EXP(M2)
         A=-PE*M1*EXPM1
         B=-PE*M2*EXPM2
         C1=-B/(A-B)
         C2= A/(A-B)
C...
C...  EQUATION (17)
         UA2=C1*EXP(M1*Z)+C2*EXP(M2*Z)
         RETURN
      END IF
      END
      FUNCTION UA3(PE,DA,UL,Z)
C...
C...  FUNCTION UA3 COMPUTES THE STEADY STATE SOLUTION FOR CASE 3
C...
      REAL M1, M2
C...
C...  EQUATION (27)
      IF(DA.EQ.0.)THEN
         UA3=1.0
         RETURN
C...  
C...  EQUATION (16)
      ELSE
         M1=(PE+SQRT(PE**2+4.0*DA))/2.0
         M2=(PE-SQRT(PE**2+4.0*DA))/2.0
C...
C...  EQUATIONS (20) AND (21) (NOTE - FOR LARGE PE, E.G., PE = 100,
C...  THESE EXPONENTIALS OVERFLOW, BUT GENERALLY, THE NUMERICAL SOLU-
C...  TION CAN BE COMPUTED)
         EXPM1=EXP(M1)
         EXPM2=EXP(M2)
         A=-(PE*M1+DA)*EXPM1
         B=-(PE*M2+DA)*EXPM2
         C1=-B/(A-B)
         C2= A/(A-B)
C...
C...  EQUATION (17)
         UA3=C1*EXP(M1*Z)+C2*EXP(M2*Z)
         RETURN
      END IF
      END
