*DECK PRO6P1
      SUBROUTINE INITAL
C...
C...  ONE DIMENSIONAL BURGERS' EQUATION WITH DIRICHLET AND NEUMANN
C...  BOUNDARY CONDITIONS, ILLUSTRATING THE NUMERICAL METHOD OF LINES
C...  INTEGRATION OF ONE DIMENSIONAL ELLIPTIC, HYPERBOLIC AND PARABOLIC
C...  PARTIAL DIFFERENTIAL EQUATIONS (PDES) - EXPLICIT TIME INTEGRATION
C...
C...  THIS PROBLEM ILLUSTRATES THE METHOD OF LINES PROGRAMMING OF A
C...  ONE DIMENSIONAL HYPERBOLIC PARABOLIC (CONVECTIVE DIFFUSION)
C...  PARTIAL DIFFERENTIAL EQUATION (PDE), THE WELL KNOWN BURGERS'
C...  EQUATION
C...
C...     U  = -U*U  + VIS*U                                     (1)
C...      T       X        XX
C...
C...  WHERE
C...
C...     U       DEPENDENT VARIABLE TO BE COMPUTED
C...
C...     T       INITIAL VALUE INDEPENDENT VARIABLE (TYPICALLY TIME)
C...
C...     X       BOUNDARY VALUE INDEPENDENT VARIABLE (TYPICALLY SPACE)
C...
C...     VIS     COEFFICIENT, TYPICALLY A VISCOSITY, WHICH CHANGES THE
C...             CHARACTER OF EQUATION (1)
C...
C...                VIS LARGE - EQUATION (1) IS STRONGLY PARABOLIC
C...
C...                VIS SMALL - EQUATION (1) IS STRONGLY HYPERBOLIC
C...
C...  A LETTER SUBSCRIPT DENOTES A PARTIAL DERIVATIVE, E.G., U  IS THE
C...                                                          T
C...  FIRST ORDER PARTIAL DERIVATIVE OF U WITH RESPECT TO T.
C...
C...  EQUATION (1) IS A WELL ESTABLISHED TEST PROBLEM FOR PDE SOFTWARE
C...  BECAUSE IT
C...
C...     (1)  HAS AN EXACT SOLUTION (A), E.G.,
C...
C...             U(X,T) = 1/(1 + EXP(X/(2*VIS) - T/(4*VIS)))    (2)
C...
C...           (A)  BENTON, E. R., AND G. W. PLATZMAN, A TABLE OF
C...           SOLUTIONS OF THE ONE DIMENSIONAL BURGERS EQUATION,
C...           QUARTERLY OF APPLIED MATHEMATICS, JULY, 1972, PP
C...           195-212.
C...
C...     (2)  IS A STRINGENT TEST PROBLEM, PARTICULARLY FOR SMALL VIS,
C...          FOR WHICH THE SOLUTION EXHIBITS STEEP, MOVING FRONTS.
C...
C...  IN THE SUBSEQUENT PROGRAMMING OF EQUATION (1), THE INITIAL AND
C...  BOUNDARY CONDITIONS HAVE BEEN CHOSEN TO BE CONSISTENT WITH
C...  EQUATION (2).  ALSO, VIS = 1.0, WHICH GIVES EQUATION (1) A
C...  MIXED HYPERBOLIC PARABOLIC CHARACTER, RESULTS IN A VERY SMOOTH
C...  SOLUTION FROM EQUATION (2), I.E., THIS CHOICE OF VIS WAS MADE
C...  TO FACILITATE THE CALCULATION OF THE NUMERICAL SOLUTION IN ORDER
C...  TO DEMONSTRATE VARIOUS COMPUTATIONAL AND PROGRAMMING PROCEDURES
C...  WITH A MODEST COMPUTING EFFORT.
C...
C...  EQUATION (1) IS INTEGRATED BY THE NUMERICAL METHOD OF LINES IN
C...  THE SUBSEQUENT PROGRAMMING FOR TWO TYPES OF BOUNDARY CONDITIONS
C...  AT X = 0 AND X = XL (IN THE SUBSEQUENT DISCUSSION, XL HAS THE
C...  NUMERICAL VALUE XL = 1)
C...
C...     (1)  DIRICHLET BOUNDARY CONDITIONS IN WHICH U(0,T) AND
C...          U(XL,T) ARE SPECIFIED
C...
C...     (2)  NEUMANN BOUNDARY CONDITIONS IN WHICH U (0,T) AND
C...          U (XL,T) ARE SPECIFIED                X
C...           X
C...
C...  ALSO, TWO APPROACHES TO THE CALCULATION OF THE DERIVATIVE U
C...                                                             XX
C...  (WHICH WILL SUBSEQUENTLY BE DESIGNATED WITH THE IN LINE NOTATION
C...  UXX) IN EQUATION (1) ARE USED
C...
C...     (1)  STAGEWISE DIFFERENTIATION, I.E., UX IS CALCULATED FROM U
C...          AND THEN UXX IS CALCULATED FROM UX
C...
C...     (2)  DIRECT DIFFERENTIATION, I.E., UXX IS CALCULATED DIRECTLY
C...          FROM U
C...
C...  THUS, THERE ARE A TOTAL OF FOUR COMBINATIONS (2 TYPES OF BOUNDARY
C...  CONDITIONS X 2 TYPES OF UXX CALCULATIONS) IN SUBROUTINES DERV1,
C...  DERV2, DERV3 AND DERV4.  COMMENTS ARE GIVEN IN EACH OF THESE
C...  SUBROUTINES EXPLAINING THE PROGRAMMING IN DETAIL.
C...
C...  SINCE THE EXACT SOLUTION, EQUATION (2), IS AVAILABLE, THE ERROR
C...  IN THE NUMERICAL SOLUTION CAN BE COMPUTED.  IN PARTICULAR, THE
C...  ERROR IN U(X=XL/2,T=1.0) IS COMPUTED FOR 11, 21, 31, 41 AND 51
C...  GRID POINTS IN X IN THE NUMERICAL METHOD OF LINES SOLUTION.  THIS
C...  ERROR IS THEN PLOTTED VS THE GRID SPACING LOGARITHMICALLY AT THE
C...  END OF THE 20TH RUN (4 COMBINATIONS MENTIONED PREVIOUSLY X 5
C...  NUMBERS OF GRID POINTS FOR EACH COMBINATION).  THESE CALCULATIONS,
C...  AND THE ASSOCIATED PLOTTING, ARE EXPLAINED IN THE COMMENTS IN
C...  SUBROUTINE PRINT.
C...
C...  THE EXACT SOLUTION IN THIS CASE IS (WITH X = XL/2 = 0.5, T = 1)
C...
C...     U(0.5,1.0) = 1/(1 + EXP(0.5/(2*VIS) - 1/(4*VIS)))
C...
C...                = 1/(1 + EXP(0)) = 1/2 = 0.5
C...
C...  NOTE THAT THIS RESULT IS INDEPENDENT OF VIS.  THUS, VIS CAN BE
C...  DECREASED TO MAKE THE NUMERICAL SOLUTION MORE DIFFICULT TO
C...  COMPUTE, I.E., THE PROBLEM IS MORE HYPERBOLIC AND LESS PARABOLIC,
C...  YET THE NUMERICAL SOLUTION SHOULD REMAIN CLOSE TO 0.5 IF REASON-
C...  ABLE ACCURACY IS TO BE MAINTAINED.
C...
C...  IN ADDITION TO THE COMMENTS IN THE FOLLOWING CODING, THE READER
C...  CAN ALSO REFER TO THE COMMENTS IN THE SPATIAL DIFFERENTIATION
C...  ROUTINES, DSS004, DSS020 AND DSS044.  IN PARTICULAR, CENTERED
C...  APPROXIMATIONS IN SUBROUTINES DSS004 AND DSS044 ARE USED TO
C...  CALCULATE UXX IN EQUATION (1), I.E., CENTERED APPROXIMATIONS
C...  SHOULD BE USED IN PARABOLIC PROBLEMS, AND NONCENTERED (FIVE-
C...  POINT BIASED UPWIND) APPROXIMATIONS IN SUBROUTINE DSS020 ARE
C...  USED TO CALCULATE UX IN EQUATION (1), I.E., NONCENTERED
C...  APPROXIMATIONS SHOULD BE USED IN HYPERBOLIC PROBLEMS.  ALL
C...  OF THESE APPROXIMATIONS ARE DERIVED FROM THE TAYLOR SERIES
C...  AS EXPLAINED IN THE THREE SUBROUTINES.
C...
C...  TWO SPECIAL CASES OF THIS PROBLEM CAN BE CONSIDERED
C...
C...     (1)  VIS = 0
C...
C...     IN THIS CASE, EQUATION (1) REDUCES TO
C...
C...     U  = -U*U  = -(U**2) /2                                (3)
C...      T       X          X
C...
C...     EQUATION (3) HAS THE ANALYTICAL SOLUTION
C...
C...     U(X,T) = X/(1 + T)                                     (4)
C...
C...     THE USE OF EQUATIONS (3) AND (4) IN THE FOLLOWING PROGRAMMING
C...     REQUIRES SOME MODIFICATION OF THE PROGRAMMING.
C...
C...     (2)  BOUNDARY VALUE ORDINARY DIFFERENTIAL EQUATION (ODE)
C...
C...     THE STEADY STATE FORM OF EQUATION (1) IS THE BOUNDARY VALUE
C...     ORDINARY DIFFERENTIAL EQUATION
C...
C...     U   - U*U  = 0                                         (5)
C...      XX      X
C...
C...     (NOTE THAT U  = 0 IN EQUATION (1)), WHICH CAN ALSO BE CONSIDER-
C...                 T
C...     AS A ONE DIMENSIONAL ELLIPTIC PDE (WITH ONE INDEPENDENT VARI-
C...     ABLE, X).
C...
C...     EQUATION (5) CAN BE SOLVED SUBJECT TO TWO BOUNDARY CONDITIONS
C...     BY ALLOWING THE SOLUTION OF EQUATION (1) TO PROCEED TO U  = 0.
C...                                                             T
C...     (AN EXACT SOLUTION CAN ALSO BE EASILY DERIVED).  IF THE BOUND-
C...     ARY CONDITIONS ARE SELECTED AS U(0,INF) = U(XL,INF) = 1, THE
C...     EXACT SOLUTION IS U(X,INF) = 1.  THIS SOLUTION CAN BE DEMON-
C...     STRATED BY ALLOWING THE NUMERICAL SOLUTION TO EQUATION (1) TO
C...     PROCEED TO T = INF, I.E., U  = 0.
C...                                T
C...
C...     THIS DISCUSSION SUGGESTS A GENERAL METHOD FOR SOLVING BOUNDARY
C...     VALUE ODES, I.E., A PSEUDO TIME DERIVATIVE IS ADDED TO THE ODES
C...     AND THE RESULTING PDES ARE INTEGRATED TO STEADY STATE FOR WHICH
C...     THE TIME DERIVATIVES ESSENTIALLY GO TO ZERO.  THIS APPROACH HAS
C...     BEEN USED EFFECTIVELY IN THE SOLUTION OF SOME RELATIVELY COM-
C...     PLEX BOUNDARY VALUE ODES (WHICH, AGAIN, CAN ALSO BE CONSIDERED
C...     ONE DIMENSIONAL ELLIPTIC PDES).  THIS APPROACH IS ALSO TERMED
C...     THE METHOD OF FALSE TRANSIENTS.
C...
C...  FINALLY, THE VERIFICATION OF EQUATION (2) AS A SOLUTION TO EQUA-
C...  TION (1) FOLLOWS:
C...
C...  LET EX = EXP(X/(2*VIS) - T/(4*VIS)).  THEN EQUATION (2) CAN BE
C...  WRITTEN AS:
C...
C...     U(X,T) = 1/(1 + EX)
C...
C...  THE VARIOUS TERMS IN EQUATION (1) ARE:
C...
C...     U  = (-1/(1 + EX)**2)*EX*(-1/(4*VIS))
C...      T
C...
C...     U*U  = (1/(1 + EX))*(-1/(1 + EX)**2)*EX*(1/(2*VIS))
C...        X
C...
C...     -VIS*U   = -VIS*((-1/(1 + EX)**2)*EX*(1/(2*VIS))**2
C...           XX
C...
C...              + EX*(1/(2*VIS)*2/(1 + EX)**3)*EX*(1/(2*VIS)))
C...
C...  SUMMATION OF THESE TERMS GIVES:
C...
C...     0 = (1/(1 + EX)**2)*EX*(1/(4*VIS))*(1 + 1)
C...
C...       - (1/(1 + EX)**3)*EX*(1/(2*VIS))*(1 + EX)
C...
C...       = 0
C...
C...  WE NOW PROCEED TO THE NUMERICAL METHOD OF LINES SOLUTION OF EQUA-
C...  TION (1) SUBJECT TO DIRICHLET AND NEUMANN BOUNDARY CONDITIONS, AND
C...  AN INITIAL CONDITION, WHICH ARE CONSISTENT WITH THE EXACT SOLU-
C...  TION, EQUATION (2).
C...
C...  ODE COMMON
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...  EQUATION PARAMETERS
      VIS=1.0
      XL =1.0E+00
C...
C...  NUMBER OF GRID POINTS
C...     DIRICHLET BOUNDARY CONDITIONS
         IF((NORUN- 1)*(NORUN- 6).EQ.0)N=11
         IF((NORUN- 2)*(NORUN- 7).EQ.0)N=21
         IF((NORUN- 3)*(NORUN- 8).EQ.0)N=31
         IF((NORUN- 4)*(NORUN- 9).EQ.0)N=41
         IF((NORUN- 5)*(NORUN-10).EQ.0)N=51
C...
C...     NEUMANN BOUNDARY CONDITIONS
         IF((NORUN-11)*(NORUN-16).EQ.0)N=11
         IF((NORUN-12)*(NORUN-17).EQ.0)N=21
         IF((NORUN-13)*(NORUN-18).EQ.0)N=31
         IF((NORUN-14)*(NORUN-19).EQ.0)N=41
         IF((NORUN-15)*(NORUN-20).EQ.0)N=51
C...
C...  INITIAL CONDITIONS (T = 0 IN EQUATION (2))
      DO 1 I=1,N
      X(I)=XL*FLOAT(I-1)/FLOAT(N-1)
      U(I)=UA(X(I),T,VIS)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...  DIRICHLET BOUNDARY CONDITIONS, STAGEWISE DIFFERENTIATION
      IF((NORUN.GE. 1).AND.(NORUN.LE. 5))CALL DERV1
C...
C...  DIRICHLET BOUNDARY CONDITIONS, DIRECT DIFFERENTIATION
      IF((NORUN.GT. 6).AND.(NORUN.LE.10))CALL DERV2
C...
C...  NEUMANN BOUNDARY CONDITIONS, STAGEWISE DIFFERENTIATION
      IF((NORUN.GE.11).AND.(NORUN.LE.15))CALL DERV3
C...
C...  NEUMANN BOUNDARY CONDITIONS, DIRECT DIFFERENTIATION
      IF((NORUN.GE.16).AND.(NORUN.LE.20))CALL DERV4
      RETURN
      END
      SUBROUTINE DERV1
C...
C...  DIRICHLET BOUNDARY CONDITIONS, STAGEWISE DIFFERENTIATION
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         U(1)=UA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         U(N)=UA(X(N),T,VIS)
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX IN UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS004(0.0E+00,XL,N,U,UX)
C...
C...     UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS004(0.0E+00,XL,N,UX,UXX)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS020(0.0E+00,XL,N,U,UX,1.0E+00)
C...
C...  BURGERS' EQUATION
C...
C...     X = 0
         UT(1)=0.0E+00
C...
C...     0 LT X LT XL
         NM1=N-1
         DO 1 I=2,NM1
         UT(I)=-U(I)*UX(I)+VIS*UXX(I)
1        CONTINUE
C...
C...     X = XL
         UT(N)=0.0E+00
      RETURN
      END
      SUBROUTINE DERV2
C...
C...  DIRICHLET BOUNDARY CONDITIONS, DIRECT DIFFERENTIATION
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         U(1)=UA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         U(N)=UA(X(N),T,VIS)
C...
C...  SPATIAL DERIVATIVES
C...
C...     UXX (BY SECOND DERIVATIVE APPROXIMATIONS)
         CALL DSS044(0.0E+00,XL,N,U,UX,UXX,1,1)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS020(0.0E+00,XL,N,U,UX,1.0E+00)
C...
C...  BURGERS' EQUATION
C...
C...     X = 0
         UT(1)=0.0E+00
C...
C...     0 LT X LT XL
         NM1=N-1
         DO 1 I=2,NM1
         UT(I)=-U(I)*UX(I)+VIS*UXX(I)
1        CONTINUE
C...
C...     X = XL
         UT(N)=0.0E+00
      RETURN
      END
      SUBROUTINE DERV3
C...
C...  NEUMANN BOUNDARY CONDITIONS, STAGEWISE DIFFERENTIATION
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX IN UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS004(0.0E+00,XL,N,U,UX)
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         UX(1)=DUA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         UX(N)=DUA(X(N),T,VIS)
C...
C...     UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS004(0.0E+00,XL,N,UX,UXX)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS020(0.0E+00,XL,N,U,UX,1.0E+00)
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         UX(1)=DUA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         UX(N)=DUA(X(N),T,VIS)
C...
C...  BURGERS' EQUATION
C...
C...     0 LE X LE XL
         DO 1 I=1,N
         UT(I)=-U(I)*UX(I)+VIS*UXX(I)
1        CONTINUE
      RETURN
      END
      SUBROUTINE DERV4
C...
C...  NEUMANN BOUNDARY CONDITIONS, DIRECT DIFFERENTIATION
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         UX(1)=DUA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         UX(N)=DUA(X(N),T,VIS)
C...
C...  SPATIAL DERIVATIVES
C...
C...     UXX (BY SECOND DERIVATIVE APPROXIMATIONS)
         CALL DSS044(0.0E+00,XL,N,U,UX,UXX,2,2)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS020(0.0E+00,XL,N,U,UX,1.0E+00)
C...
C...     X = 0 (X = 0 IN EQUATION (2))
         UX(1)=DUA(X(1),T,VIS)
C...
C...     X = XL (X = XL IN EQUATION (2))
         UX(N)=DUA(X(N),T,VIS)
C...
C...  BURGERS' EQUATION
C...
C...     0 LE X LE XL
         DO 1 I=1,N
         UT(I)=-U(I)*UX(I)+VIS*UXX(I)
1        CONTINUE
      RETURN
      END
      FUNCTION UA(X,T,VIS)
C...
C...  FUNCTION UA IMLEMENTS THE ANALYTICAL SOLUTION, EQUATION (2)
C...
      EX=EXP(X/(2.0E+00*VIS)-T/(4.0E+00*VIS))
      UA=1.0E+00/(1.0E+00+EX)
      RETURN
      END
      FUNCTION DUA(X,T,VIS)
C...
C...  FUNCTION DUA IMPLEMENTS THE DERIVATIVE OF THE ANALYTICAL SOLUTION,
C...  EQUATION (2), WITH RESPECT TO X
C...
      EX=EXP(X/(2.0E+00*VIS)-T/(4.0E+00*VIS))
      DUA=-1.0E+00*(1.0E+00+EX)**(-2)*EX*(1.0E+00/(2.0E+00*VIS))
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     U(51)
     2      /F/    UT(51)
     3      /S/    UX(51),   UXX(51)
     4      /R/       VIS,        XL,     X(51)
     5      /I/         N,        IP
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION DXP(5), EP1(2,5), EP2(2,5)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(11X,'T',7X,'U(XL/2,T)',6X,'U(XL/2,T)',8X,'ERROR',/,
     1              19X,'  (NUM)  ',6X,'  (ANAL) ')
C...
C...  COMPUTE THE ANALYTICAL SOLUTION AND THE ERROR IN THE NUMERICAL
C...  SOLUTION AT X = XL/2
      N2=N/2+1
      XL2=X(N2)
      UANAL=UA(XL2,T,VIS)
      ERROR=U(N2)-UANAL
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS, AND THE ERROR
      WRITE(NO,2)T,U(N2),UANAL,ERROR
2     FORMAT(F12.2,2F15.6,E16.4)
C...
C...  STORE THE SOLUTION U(X=0.5,T=1.0)FOR PLOTTING
C...
C...     THE ERRORS FOR NORUN = 1 TO 5 (NPLOT = 1), AND NORUN = 6 TO
C...     10 (NPLOT = 2) ARE PLOTTED VS THE GRID SPACING ON THE SAME
C...     LOG-LOG PLOT (DIRICHLET BOUNDARY CONDITIONS)
         IF((T.GT.0.95E+00).AND.(T.LT.1.05E+00))THEN
C...
C...        NPLOT=1
            IF((NORUN.GE. 1).AND.(NORUN.LE. 5))THEN
               NPLOT=1
               DXP(NORUN)=ALOG10(XL/FLOAT(N-1))
               EP1(NPLOT,NORUN-0)=ALOG10(ABS(ERROR))
C...
C...        NPLOT=2
            ELSE IF((NORUN.GE. 6).AND.(NORUN.LE.10))THEN
               NPLOT=2
               EP1(NPLOT,NORUN-5)=ALOG10(ABS(ERROR))
            END IF
C...
C...     THE ERRORS FOR NORUN = 11 TO 15 (NPLOT = 1), AND NORUN = 16 TO
C...     10 (NPLOT = 2) ARE PLOTTED VS THE GRID SPACING ON THE SAME
C...     LOG-LOG PLOT (NEUMANN BOUNDARY CONDITIONS)
C...
C...        NPLOT = 1
            IF((NORUN.GE.11).AND.(NORUN.LE.15))THEN
               NPLOT=1
               EP2(NPLOT,NORUN-10)=ALOG10(ABS(ERROR))
C...
C...        NPLOT=2
            ELSE IF((NORUN.GE.16).AND.(NORUN.LE.20))THEN
               NPLOT=2
               EP2(NPLOT,NORUN-15)=ALOG10(ABS(ERROR))
            END IF
         END IF
C...
C...  PLOT THE ERROR VS GRID SPACING AT THE END OF FOUR SERIES OF FIVE
C...  RUNS
      IF((T.GT.1.95E+00).AND.(T.LT.2.05E+00).AND.(NORUN.EQ.20))THEN
C...
C...     DIRICHLET BC
         CALL TPLOTS(2,5,DXP,EP1)
C...
C...     PRINT A LABEL FOR THE PLOT
         WRITE(NO,11)
11       FORMAT(
     1 ' DIRICHLET BC, LOG(ERROR) VS LOG(DX)')
C...
C...     NEUMANN BC
         CALL TPLOTS(2,5,DXP,EP2)
C...
C...     PRINT A LABEL FOR THE PLOT
         WRITE(NO,12)
12       FORMAT(
     1 ' NEUMANN BC, LOG(ERROR) VS LOG(DX)')
      END IF
      RETURN
      END
