*DECK PRO6P3
      SUBROUTINE INITAL
C...
C...  THREE DIMENSIONAL BURGERS' EQUATION WITH DIRICHLET AND NEUMANN
C...  BOUNDARY CONDITIONS, ILLUSTRATING THE NUMERICAL METHOD OF LINES
C...  INTEGRATION OF THREE DIMENSIONAL ELLIPTIC, HYPERBOLIC AND PARA-
C...  BOLIC PARTIAL DIFFERENTIAL EQUATIONS (PDES)
C...
C...  THIS PROBLEM ILLUSTRATES THE METHOD OF LINES PROGRAMMING OF A
C...  THREE DIMENSIONAL HYPERBOLIC PARABOLIC (CONVECTIVE DIFFUSION)
C...  PARTIAL DIFFERENTIAL EQUATION (PDE), THE WELL KNOWN BURGERS'
C...  EQUATION
C...
C...     U  = -U*U  + VIS*U   - U*U  + VIS*U                   (1)
C...      T       X        XX      Y        YY
C...
C...  WHERE
C...
C...     U       DEPENDENT VARIABLE TO BE COMPUTED
C...
C...     T       INITIAL VALUE INDEPENDENT VARIABLE (TYPICALLY TIME)
C...
C...     X,Y,Z   BOUNDARY VALUE INDEPENDENT VARIABLES (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, E.G.,
C...
C...          U(X,Y,T) =
C...
C...          1/(1 + EXP(X/(2*VIS) + Y/(2*VIS) + Z/(2*VIS)      (2)
C...
C...               - 3*T/(4*VIS)))
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 (1) X = AND X = XL, (2) Y = 0 AND Y = YL AND (3) Z = 0 AND
C...  Z = ZL (IN THE SUBSEQUENT DISCUSSION, XL, YL AND ZL HAVE THE
C...  NUMERICAL VALUES XL = YL = ZL = 1)
C...
C...     (1)  DIRICHLET BOUNDARY CONDITIONS IN WHICH U(0,Y,Z,T),
C...          U(XL,Y,Z,T), U(X,0,Z,T), U(X,YL,Z,T), U(X,Y,0,T) AND
C...          U(X,Y,ZL,T) ARE SPECIFIED
C...
C...     (2)  NEUMANN BOUNDARY CONDITIONS IN WHICH U (0,Y,Z,T),
C...                                                X
C...          U (XL,Y,Z,T), U (X,0,Z,T), U (X,YL,Z,T), U (X,Y,0,T)
C...           X             Y            Y             Z
C...          AND U (X,Y,ZL,T) ARE SPECIFIED
C...               Z
C...
C...  IN THE SUBSEQUENT DISCUSSION, THE DERIVATIVES IN EQUATION (1)
C...  WILL BE DESIGNATED WITH IN-LINE NOTATION, E.G., U   WILL BE
C...                                                   XX
C...  DESIGNATED UXX.
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,Y=YL/2,Z=ZL/2,T=1) IS COMPUTED FROM THE FOLLOW-
C...  ING APPLICATION OF EQUATION (2)
C...
C...  U(0.5,0.5,0.5,1) = 1/(1 + EXP(0.5/(2*VIS) + 0.5/(2*VIS)
C...
C...                              + 0.5/(2*VIS) - 3*1/(4*VIS)))
C...
C...                   = 1/(1 + EXP(0)) = 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...  SUBROUTINE, DSS036.  IN PARTICULAR, CENTERED APPROXIMATIONS IN
C...  SUBROUTINE DSS036 ARE USED TO CALCULATE UXX AND UYY IN EQUATION
C...  (1), I.E., CENTERED APPROXIMATIONS SHOULD BE USED IN PARABOLIC
C...  PROBLEMS, AND NONCENTERED (FIVE POINT BIASED UPWIND) APPROXI-
C...  MATIONS IN SUBROUTINE DSS036 ARE USED TO CALCULATE UX AND UY IN
C...  EQUATION (1), I.E., NONCENTERED APPROXIMATIONS SHOULD BE USED
C...  IN HYPERBOLIC PROBLEMS.  ALL OF THESE APPROXIMATIONS ARE DERIVED
C...  FROM THE TAYLOR SERIES AS EXPLAINED IN DSS036, AND SUBORDINATE
C...  ROUTINES DSS004 AND DSS020.
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*U  - U*U
C...      T        X      Y      Z
C...                                                            (3)
C...          - (U**2) /2 - (U**2) /2 - (U**2) /2
C...                  X           Y           Z
C...
C...     EQUATION (3) HAS THE ANALYTICAL SOLUTION
C...
C...     U(X,Y,Z,T) = (X + Y + Z)/(1 + 3*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 PARTIAL DIFFERENTIAL EQUATION (PDE)
C...
C...     THE STEADY STATE FORM OF EQUATION (1) IS THE BOUNDARY VALUE
C...     PARTIAL DIFFERENTIAL EQUATION
C...
C...     U   - U*U  + U   - U*U  + U   - U*U  = 0               (5)
C...      XX      X    YY      Y    ZZ      Z
C...
C...     (NOTE THAT U  = 0 IN EQUATION (1)), WHICH CAN ALSO BE CONSIDER-
C...                 T
C...     ED AS A THREE DIMENSIONAL ELLIPTIC PDE (WITH INDEPENDENT VARI-
C...     ABLES X, Y AND Z).  THUS, THE FOLLOWING PROGRAMMING ILLUSTRATES
C...     A METHOD FOR INTEGRATING ELLIPTIC PDES AS WELL AS HYPERBOLIC
C...     AND PARABOLIC PDES (BY ALLOWING THE SOLUTION OF EQUATION (1)
C...     TO PROCEED TO U  = 0, FOR WHICH EQUATION (2), WITH LARGE T,
C...                    T
C...     REDUCES TO:
C...
C...        U(X,Y,Z,INF) = 1
C...
C...     THIS DISCUSSION SUGGESTS A GENERAL METHOD FOR SOLVING BOUNDARY
C...     VALUE PDES, I.E., A PSEUDO TIME DERIVATIVE IS ADDED TO THE PDES
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 PDES.  THIS APPROACH IS ALSO TERMED THE
C...     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) + Y/(2*VIS) + Z/(2*VIS) - 3*T/(4*VIS)).
C...  THEN EQUATION (2) CAN BE WRITTEN AS:
C...
C...     U(X,Y,Z,T) = 1/(1 + EX)
C...
C...  THE VARIOUS TERMS IN EQUATION (1) ARE:
C...
C...     U  = (-1/(1 + EX)**2)*EX*(-3/(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...     U*U  = (1/(1 + EX))*(-1/(1 + EX)**2)*EX*(1/(2*VIS))
C...        Y
C...
C...     -VIS*U   = -VIS*((-1/(1 + EX)**2)*EX*(1/(2*VIS))**2
C...           YY
C...
C...              + EX*(1/(2*VIS)*2/(1 + EX)**3)*EX*(1/(2*VIS)))
C...
C...     U*U  = (1/(1 + EX))*(-1/(1 + EX)**2)*EX*(1/(2*VIS))
C...        Z
C...
C...     -VIS*U   = -VIS*((-1/(1 + EX)**2)*EX*(1/(2*VIS))**2
C...           ZZ
C...
C...  SUMMATION OF THESE TERMS GIVES:
C...
C...     0 = (1/(1 + EX)**2)*EX*(1/(4*VIS))*(3 + 1 + 1 + 1)
C...
C...       - (1/(1 + EX)**3)*EX*(1/(2*VIS))*(1 + 1 + 1)*(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(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...  EQUATION PARAMETERS
      VIS=1.0E+00
      XL=1.0E+00
      YL=1.0E+00
      ZL=1.0E+00
C...
C...  NUMBER OF GRID POINTS
      NX=6
      NY=6
      NZ=6
C...
C...  INITIAL CONDITIONS (T = 0 IN EQUATION (2))
      DO 1 I=1,NX
      X(I)=XL*FLOAT(I-1)/FLOAT(NX-1)
      DO 1 J=1,NY
      Y(J)=YL*FLOAT(J-1)/FLOAT(NY-1)
      DO 1 K=1,NZ
      Z(K)=ZL*FLOAT(K-1)/FLOAT(NZ-1)
      U(I,J,K)=UA(X(I),Y(J),Z(K),T,VIS)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...  DIRICHLET BOUNDARY CONDITIONS WITH (1) UXX, UYY AND UZZ CALCULA-
C...  TED BY CENTERED APPROXIMATIONS AND (2) UX, UY AND UZ CALCULATED
C...  BY BIASED UPWIND APPROXIMATIONS
      IF(NORUN.EQ.1)CALL DERV1
C...
C...  DIRICHLET BOUNDARY CONDITIONS WITH UXX, UYY, UZZ, UX, UY AND UZ
C...  CALCULATED BY CENTERED APPROXIMATIONS
      IF(NORUN.EQ.2)CALL DERV2
C...
C...  NEUMANN BOUNDARY CONDITIONS WITH (1) UXX, UYY AND UZZ CALCULA-
C...  TED BY CENTERED APPROXIMATIONS AND (2) UX, UY AND UZ CALCULATED
C...  BY BIASED UPWIND APPROXIMATIONS
      IF(NORUN.EQ.3)CALL DERV3
C...
C...  NEUMANN BOUNDARY CONDITIONS WITH UXX, UYY, UZZ, UX, UY AND UZ
C...  CALCULATED BY CENTERED APPROXIMATIONS
      IF(NORUN.EQ.4)CALL DERV4
      RETURN
      END
      SUBROUTINE DERV1
C...
C...  DIRICHLET BOUNDARY CONDITIONS
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...     X = 0, 0 LE Y LE YL, 0 LE Z LE ZL (X = 0 IN EQUATION (2))
         DO 1 J=1,NY
         DO 1 K=1,NZ
         U(1,J,K)=UA(X(1),Y(J),Z(K),T,VIS)
1        CONTINUE
C...
C...     X = XL, 0 LE Y LE YL, 0 LE Z LE ZL (X = XL IN EQUATION (2))
         DO 2 J=1,NY
         DO 2 K=1,NZ
         U(NX,J,K)=UA(X(NX),Y(J),Z(K),T,VIS)
2        CONTINUE
C...
C...     Y = 0, 0 LE X LE XL, 0 LE Z LE ZL (Y = 0 IN EQUATION (2))
         DO 3 I=1,NX
         DO 3 K=1,NZ
         U(I,1,K)=UA(X(I),Y(1),Z(K),T,VIS)
3        CONTINUE
C...
C...     Y = YL, 0 LE X LE XL, 0 LE Z LE ZL (Y = YL IN EQUATION (2))
         DO 4 I=1,NX
         DO 4 K=1,NZ
         U(I,NY,K)=UA(X(I),Y(NY),Z(K),T,VIS)
4        CONTINUE
C...
C...     Z = 0, 0 LE X LE XL, 0 LE Y LE YL (Z = 0 IN EQUATION (2))
         DO 5 I=1,NX
         DO 5 J=1,NY
         U(I,J,1)=UA(X(I),Y(J),Z(1),T,VIS)
5        CONTINUE
C...
C...     Z = ZL, 0 LE X LE XL, 0 LE Y LE YL (Z = ZL EQUATION (2))
         DO 6 I=1,NX
         DO 6 J=1,NY
         U(I,J,NZ)=UA(X(I),Y(J),Z(NZ),T,VIS)
6        CONTINUE
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX IN UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,0.E+00)
C...
C...     UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,UX,UXX,0.E+00)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,1.0E+00)
C...
C...     UY IN UYY (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,0.E+00)
C...
C...     UYY (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,UY,UYY,0.E+00)
C...
C...     UY IN U*UY (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,1.0E+00)
C...
C...     UZ IN UZZ (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,0.E+00)
C...
C...     UZZ (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,UZ,UZZ,0.E+00)
C...
C...     UZ IN U*UZ (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,1.0E+00)
C...
C...  BURGERS' EQUATION
C...
         DO 10 I=1,NX
         DO 10 J=1,NY
         DO 10 K=1,NZ
         UT(I,J,K)=-U(I,J,K)*UX(I,J,K)+VIS*UXX(I,J,K)
     1             -U(I,J,K)*UY(I,J,K)+VIS*UYY(I,J,K)
     2             -U(I,J,K)*UZ(I,J,K)+VIS*UZZ(I,J,K)
10       CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
C...
C...  DIRICHLET BOUNDARY CONDITIONS
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...     X = 0, 0 LE Y LE YL, 0 LE Z LE ZL (X = 0 IN EQUATION (2))
         DO 1 J=1,NY
         DO 1 K=1,NZ
         U(1,J,K)=UA(X(1),Y(J),Z(K),T,VIS)
1        CONTINUE
C...
C...     X = XL, 0 LE Y LE YL, 0 LE Z LE ZL (X = XL IN EQUATION (2))
         DO 2 J=1,NY
         DO 2 K=1,NZ
         U(NX,J,K)=UA(X(NX),Y(J),Z(K),T,VIS)
2        CONTINUE
C...
C...     Y = 0, 0 LE X LE XL, 0 LE Z LE ZL (Y = 0 IN EQUATION (2))
         DO 3 I=1,NX
         DO 3 K=1,NZ
         U(I,1,K)=UA(X(I),Y(1),Z(K),T,VIS)
3        CONTINUE
C...
C...     Y = YL, 0 LE X LE XL, 0 LE Z LE ZL (Y = YL IN EQUATION (2))
         DO 4 I=1,NX
         DO 4 K=1,NZ
         U(I,NY,K)=UA(X(I),Y(NY),Z(K),T,VIS)
4        CONTINUE
C...
C...     Z = 0, 0 LE X LE XL, 0 LE Y LE YL (Z = 0 IN EQUATION (2))
         DO 5 I=1,NX
         DO 5 J=1,NY
         U(I,J,1)=UA(X(I),Y(J),Z(1),T,VIS)
5        CONTINUE
C...
C...     Z = ZL, 0 LE X LE XL, 0 LE Y LE YL (Z = ZL EQUATION (2))
         DO 6 I=1,NX
         DO 6 J=1,NY
         U(I,J,NZ)=UA(X(I),Y(J),Z(NZ),T,VIS)
6        CONTINUE
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,0.E+00)
C...
C...     UXX BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,UX,UXX,0.E+00)
C...
C...     UY BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,0.E+00)
C...
C...     UYY BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,UY,UYY,0.E+00)
C...
C...     UZ BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,0.E+00)
C...
C...     UZZ BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,UZ,UZZ,0.E+00)
C...
C...  BURGERS' EQUATION
C...
         DO 10 I=1,NX
         DO 10 J=1,NY
         DO 10 K=1,NZ
         UT(I,J,K)=-U(I,J,K)*UX(I,J,K)+VIS*UXX(I,J,K)
     1             -U(I,J,K)*UY(I,J,K)+VIS*UYY(I,J,K)
     2             -U(I,J,K)*UZ(I,J,K)+VIS*UZZ(I,J,K)
10       CONTINUE
      RETURN
      END
      SUBROUTINE DERV3
C...
C...  NEUMANN BOUNDARY CONDITIONS
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX IN UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,0.E+00)
C...
C...     X = 0, 0 LE Y LE YL, 0 LE Z LE ZL (X = 0 IN EQUATION (2))
         DO 1 J=1,NY
         DO 1 K=1,NZ
         UX(1,J,K)=DUA(X(1),Y(J),Z(K),T,VIS)
1        CONTINUE
C...
C...     X = XL, 0 LE Y LE YL, 0 LE Z LE ZL (X = XL IN EQUATION (2))
         DO 2 J=1,NY
         DO 2 K=1,NZ
         UX(NX,J,K)=DUA(X(NX),Y(J),Z(K),T,VIS)
2        CONTINUE
C...
C...     UXX (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,UX,UXX,0.E+00)
C...
C...     UX IN U*UX (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,1.0E+00)
C...
C...     X = 0, 0 LE Y LE YL, 0 LE Z LE ZL (X = 0 IN EQUATION (2))
         DO 7 J=1,NY
         DO 7 K=1,NZ
         UX(1,J,K)=DUA(X(1),Y(J),Z(K),T,VIS)
7        CONTINUE
C...
C...     X = XL, 0 LE Y LE YL, 0 LE Z LE ZL (X = XL IN EQUATION (2))
         DO 8 J=1,NY
         DO 8 K=1,NZ
         UX(NX,J,K)=DUA(X(NX),Y(J),Z(K),T,VIS)
8        CONTINUE
C...
C...     UY IN UYY (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,0.E+00)
C...
C...     Y = 0, 0 LE X LE XL, 0 LE Z LE ZL (Y = 0 IN EQUATION (2))
         DO 3 I=1,NX
         DO 3 K=1,NZ
         UY(I,1,K)=DUA(X(I),Y(1),Z(K),T,VIS)
3        CONTINUE
C...
C...     Y = YL, 0 LE X LE XL, 0 LE Z LE ZL (Y = YL IN EQUATION (2))
         DO 4 I=1,NX
         DO 4 K=1,NZ
         UY(I,NY,K)=DUA(X(I),Y(NY),Z(K),T,VIS)
4        CONTINUE
C...
C...     UYY (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,UY,UYY,0.E+00)
C...
C...     UY IN U*UY (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,1.0E+00)
C...
C...     Y = 0, 0 LE X LE XL, 0 LE Z LE ZL (Y = 0 IN EQUATION (2))
         DO 9 I=1,NX
         DO 9 K=1,NZ
         UY(I,1,K)=DUA(X(I),Y(1),Z(K),T,VIS)
9        CONTINUE
C...
C...     Y = YL, 0 LE X LE XL, 0 LE Z LE ZL (Y = YL IN EQUATION (2))
         DO 10 I=1,NX
         DO 10 K=1,NZ
         UY(I,NY,K)=DUA(X(I),Y(NY),Z(K),T,VIS)
10       CONTINUE
C...
C...     UZ IN UZZ (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,0.E+00)
C...
C...     Z = 0, 0 LE X LE XL, 0 LE Y LE YL (Z = 0 IN EQUATION (2))
         DO 5 I=1,NX
         DO 5 J=1,NY
         UZ(I,J,1)=DUA(X(I),Y(J),Z(1),T,VIS)
5        CONTINUE
C...
C...     Z = ZL, 0 LE X LE XL, 0 LE Y LE YL (Z = ZL EQUATION (2))
         DO 6 I=1,NX
         DO 6 J=1,NY
         UZ(I,J,NZ)=DUA(X(I),Y(J),Z(NZ),T,VIS)
6        CONTINUE
C...
C...     UZZ (BY FIVE POINT CENTERED APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,UZ,UZZ,0.E+00)
C...
C...     UZ IN U*UZ (BY FIVE POINT BIASED UPWIND APPROXIMATIONS)
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,1.0E+00)
C...
C...     Z = 0, 0 LE X LE XL, 0 LE Y LE YL (Z = 0 IN EQUATION (2))
         DO 11 I=1,NX
         DO 11 J=1,NY
         UZ(I,J,1)=DUA(X(I),Y(J),Z(1),T,VIS)
11       CONTINUE
C...
C...     Z = ZL, 0 LE X LE XL, 0 LE Y LE YL (Z = ZL EQUATION (2))
         DO 12 I=1,NX
         DO 12 J=1,NY
         UZ(I,J,NZ)=DUA(X(I),Y(J),Z(NZ),T,VIS)
12       CONTINUE
C...
C...  BURGERS' EQUATION
C...
         DO 13 I=1,NX
         DO 13 J=1,NY
         DO 13 K=1,NZ
         UT(I,J,K)=-U(I,J,K)*UX(I,J,K)+VIS*UXX(I,J,K)
     1             -U(I,J,K)*UY(I,J,K)+VIS*UYY(I,J,K)
     2             -U(I,J,K)*UZ(I,J,K)+VIS*UZZ(I,J,K)
13       CONTINUE
      RETURN
      END
      SUBROUTINE DERV4
C...
C...  NEUMANN BOUNDARY CONDITIONS
C...
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...  SPATIAL DERIVATIVES
C...
C...     UX BY FIVE POINT CENTERED APPRXOXIMATIONS
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,U,UX,0.E+00)
C...
C...     X = 0, 0 LE Y LE YL, 0 LE Z LE ZL (X = 0 IN EQUATION (2))
         DO 1 J=1,NY
         DO 1 K=1,NZ
         UX(1,J,K)=DUA(X(1),Y(J),Z(K),T,VIS)
1        CONTINUE
C...
C...     X = XL, 0 LE Y LE YL, 0 LE Z LE ZL (X = XL IN EQUATION (2))
         DO 2 J=1,NY
         DO 2 K=1,NZ
         UX(NX,J,K)=DUA(X(NX),Y(J),Z(K),T,VIS)
2        CONTINUE
C...
C...     UXX BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,XL,NX,NY,NZ,1,UX,UXX,0.E+00)
C...
C...     UY BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,U,UY,0.E+00)
C...
C...     Y = 0, 0 LE X LE XL, 0 LE Z LE ZL (Y = 0 IN EQUATION (2))
         DO 3 I=1,NX
         DO 3 K=1,NZ
         UY(I,1,K)=DUA(X(I),Y(1),Z(K),T,VIS)
3        CONTINUE
C...
C...     Y = YL, 0 LE X LE XL, 0 LE Z LE ZL (Y = YL IN EQUATION (2))
         DO 4 I=1,NX
         DO 4 K=1,NZ
         UY(I,NY,K)=DUA(X(I),Y(NY),Z(K),T,VIS)
4        CONTINUE
C...
C...     UYY BY FIVE POINT CENTERED APPROXIMATIONS
         CALL DSS036(0.E+00,YL,NX,NY,NZ,2,UY,UYY,0.E+00)
C...
C...     UZ BY FIVE POINT CENTERED APPRXIMATIONS
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,U,UZ,0.E+00)
C...
C...     Z = 0, 0 LE X LE XL, 0 LE Y LE YL (Z = 0 IN EQUATION (2))
         DO 5 I=1,NX
         DO 5 J=1,NY
         UZ(I,J,1)=DUA(X(I),Y(J),Z(1),T,VIS)
5        CONTINUE
C...
C...     Z = ZL, 0 LE X LE XL, 0 LE Y LE YL (Z = ZL EQUATION (2))
         DO 6 I=1,NX
         DO 6 J=1,NY
         UZ(I,J,NZ)=DUA(X(I),Y(J),Z(NZ),T,VIS)
6        CONTINUE
C...
C...     UZZ BY FIVE POINT APPROXIMATIONS
         CALL DSS036(0.E+00,ZL,NX,NY,NZ,3,UZ,UZZ,0.E+00)
C...
C...  BURGERS' EQUATION
C...
         DO 10 I=1,NX
         DO 10 J=1,NY
         DO 10 K=1,NZ
         UT(I,J,K)=-U(I,J,K)*UX(I,J,K)+VIS*UXX(I,J,K)
     1             -U(I,J,K)*UY(I,J,K)+VIS*UYY(I,J,K)
     2             -U(I,J,K)*UZ(I,J,K)+VIS*UZZ(I,J,K)
10       CONTINUE
      RETURN
      END
      FUNCTION UA(X,Y,Z,T,VIS)
C...
C...  FUNCTION UA IMPLEMENTS THE ANALYTICAL SOLUTION, EQUATION (2)
C...
      EX=EXP(X/(2.0E+00*VIS)
     1      +Y/(2.0E+00*VIS)
     2      +Z/(2.0E+00*VIS)
     3      -3.0E+00*T/(4.0E+00*VIS))
      UA=1.0E+00/(1.0E+00+EX)
      RETURN
      END
      FUNCTION DUA(X,Y,Z,T,VIS)
C...
C...  FUNCTION DUA IMPLEMENTS THE DERIVATIVE OF THE ANALYTICAL SOLUTION,
C...  EQUATION (2), WITH RESPECT TO X, Y AND Z
C...
      EX=EXP(X/(2.0E+00*VIS)
     1      +Y/(2.0E+00*VIS)
     2      +Z/(2.0E+00*VIS)
     3      -3.0E+00*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(6,6,6)
     2      /F/ UT(6,6,6)
     3      /S/ UX(6,6,6), UY(6,6,6), UZ(6,6,6),
     4         UXX(6,6,6),UYY(6,6,6),UZZ(6,6,6)
     5      /R/        XL,        YL,        ZL,
     6               X(6),      Y(6),      Z(6),
     7                VIS
     8      /I/        NX,        NY,        NZ,
     9                 IP
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(11X,'T',9X,'U(0.6,0.6,0.6,T)',
     1               5X,'U(0.6,0.6,0.6,T)',5X,'ERROR',/,
     2              25X,'  (NUM)  ',12X,'  (ANAL) ')
C...
C...  COMPUTE THE ANALYTICAL SOLUTION AND THE ERROR IN THE NUMERICAL
C...  SOLUTION.  NOTE THAT FOR AN EVEN NUMBER OF POINTS, E.G., 6, THE
C...  MID POINTS X = XL/2, Y = YL/2, Z = ZL/2 ARE NOT AVAILABLE.
C...  THEREFORE, THE POINTS X = 0.6*XL, Y = 0.6*YL, Z = 0.6*ZL ARE
C...  USED FOR WHICH THE EXACT SOLUTION WITH T = 1.2 IS:
C...
C...  U(0.6,0.6,0.6,1.2) =
C...
C...     = 1/(1 + EXP(0.6/(2*VIS) + 0.6/(2*VIS) + 0.6/(2*VIS)
C...
C...     - 3*1.2/(4*VIS)) = 1/(1 + EXP(0)) = 0.5
C...
      NX2=NX/2+1
      NY2=NY/2+1
      NZ2=NZ/2+1
      XL2=X(NX2)
      YL2=Y(NY2)
      ZL2=Z(NZ2)
      UANAL=UA(XL2,YL2,ZL2,T,VIS)
      ERROR=U(NX2,NY2,NZ2)-UANAL
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS, AND THE ERROR
      WRITE(NO,2)T,U(NX2,NY2,NZ2),UANAL,ERROR
2     FORMAT(F12.2,2F21.6,E16.4)
      RETURN
      END
