*DECK DSS002
      SUBROUTINE DSS002(XL,XU,N,U,UX)
C...
C...  SUBROUTINE DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU
C...
C...  ARGUMENT LIST
C...
C...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
C...
C...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
C...
C...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
C...             BOUNDARY POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
C...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
C...             TO BE COMPUTED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
C...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
C...             (OUTPUT)
C...
C...  SUBROUTINE DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM THE
C...  CLASSICAL THREE-POINT, SECOND-ORDER FINITE DIFFERENCE APPROXI-
C...  TIONS
C...
C...                                       2
C...  U1  = (1/2DX)(-3U1 + 4U2 - U3) + O(DX ) (LEFT BOUNDARY,     (1)
C...    X                                         X = XL)
C...
C...                                   2
C...  UI  = (1/2DX)(UI+1 - UI-1) + O(DX ) (INTERIOR POINT,        (2)
C...    X                                   X NE XL, XU)
C...
C...                                          2
C...  UN  = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX ) (RIGHT BOUNDARY, (3)
C...    X                                            X = XU)
C...
C...  EQUATIONS (1) TO (3) APPLY OVER A GRID IN X WITH CORRESPONDING
C...  VALUES OF THE FUNCTION U(X) REPRESENTED AS
C...
C...   U1      U2       U3         UI        UN-2      UN-1    UN
C...
C...  X=XL  X=XL+DX  X=XL+2DX ... X=XI ... X=XU-2DX  X=XU-DX  X=XU
C...
C...  THE ORIGIN OF EQUATIONS (1) TO (3) IS OUTLINED BELOW.
C...
C...  CONSIDER THE FOLLOWING POLYNOMIAL IN X OF ARBITRARY ORDER
C...
C...                                     2             3
C...  U(X) = A0 + A1(X - X0) + A2(X - X0)  + A3(X - X0)  + ....   (4)
C...
C...  WE SEEK THE VALUES OF THE COEFFICIENTS A0, A1, A2, ... FOR A
C...  PARTICULAR FUNCTION U(X).  IF X = X0 IS SUBSTITUTED IN EQUATION
C...  (4), WE HAVE IMMEDIATELY A0 = U(X0).  NEXT, IF EQUATION (4) IS
C...  DIFFERENTIATED WITH RESPECT TO X,
C...
C...                                                   2
C...  DU(X)/DX = U (X) = A1 + 2A2(X - X0) + 3A3(X - X0)  + ...    (5)
C...              X
C...
C...  AGAIN, WITH X = X0, A1 = DU(X0)/DX = U (X0).  DIFFERENTIATION
C...                                        X
C...  OF EQUATION (5) IN TURN GIVES
C...
C...  D2U(X)/DX2 = U  (X) = 2A2 + 6A3(X - X0) + ...
C...                2X
C...
C...  AND FOR X = X0, A2 = U  (X0)/2F (2F = 1*2, I.E., 2 FACTORIAL).
C...                        2X
C...
C...  WE CAN CONTINUE THIS PROCESS OF DIFFERENTIATION FOLLOWED BY THE
C...  SUBSTITUTION X = X0 TO OBTAIN THE SUCCESSIVE COEFFICIENTS IN
C...  EQUATION (4), A3, A4, ...  FINALLY, SUBSTITUTION OF THESE CO-
C...  EFFICIENTS IN EQUATION (4) GIVES
C...
C...                                                 2
C...  U(X) = U(X0) + U (X0)(X - X0) + U  (X0)(X - X0)  +
C...                  X       1F       2X       2F
C...                                                              (6)
C...                                3                  4
C...                 U  (X0)(X - X0)  + U  (X0)(X - X0)  + ...
C...                  3X       3F        4X       4F
C...
C...  THE CORRESPONDENCE BETWEEN EQUATION (6) AND THE WELL-KNOWN
C...  TAYLOR SERIES SHOULD BE CLEAR.  THUS THE EXPANSION OF A
C...  FUNCTION, U(X), AROUND A NEIGHBORING POINT X0 IN TERMS OF U(X0)
C...  AND THE DERIVATIVES OF U(X) AT X = X0 IS EQUIVALENT TO APPROXI-
C...  MATING U(X) NEAR X0 BY A POLYNOMIAL.
C...
C...  EQUATION (6) IS THE STARTING POINT FOR THE DERIVATION OF THE
C...  CLASSICAL FINITE DIFFERENCE APPROXIMATIONS OF DERIVATIVES SUCH
C...  AS THE THREE-POINT FORMULAS OF EQUATIONS (1), (2) AND (3).  WE
C...  WILL NOW CONSIDER THE DERIVATION OF THESE THREE-POINT FORMULAS
C...  IN A STANDARD FORMAT WHICH CAN THEN BE EXTENDED TO HIGHER
C...  MULTI-POINT FORMULAS IN OTHER SUBROUTINES, E.G., FIVE-POINT
C...  FORMULAS IN SUBROUTINE DSS004.
C...
C...  THREE-POINT FORMULAS
C...
C...     (1)  LEFT END, POINT I = 1
C...
C...  IF EQUATION (6) IS WRITTEN AROUND THE POINT X = XL FOR X = XL +
C...  DX AND X = XL + 2DX, FOR WHICH THE CORRESPONDING VALUES OF U(X)
C...  ARE U1, U2 AND U3 (U1 AND U2 ARE SEPARATED WITH RESPECT TO X BY
C...  DISTANCE DX AS ARE U2 AND U3, I.E., WE ASSUME A UNIFORM GRID
C...  SPACING, DX, FOR INDEPENDENT VARIABLE X)
C...
C...                                2            3
C...  U2 = U1 + U1 ( DX) + U1  ( DX)  + U1  ( DX)  + ...          (7)
C...              X  1F      2X  2F       3X  3F
C...
C...                                2            3
C...  U3 = U1 + U1 (2DX) + U1  (2DX)  + U1  (2DX)  + ...          (8)
C...              X  1F      2X  2F       3X  3F
C...
C...  WE CAN NOW TAKE A LINEAR COMBINATION OF EQUATIONS (7) AND (8)
C...  BY FIRST MULTIPLYING EQUATION (7) BY A CONSTANT, A, AND EQUA-
C...  TION (8) BY CONSTANT B
C...
C...                                  2           3
C...  A(U2 = U1 + U1 ( DX) + U1  ( DX) + U1  ( DX) + ...)         (9)
C...                X  1F      2X  2F      3X  3F
C...
C...                                  2           3
C...  B(U3 = U1 + U1 (2DX) + U1  (2DX) + U1  (2DX) + ...)        (10)
C...                X  1F      2X  2F      3X  3F
C...
C...  CONSTANTS A AND B ARE THEN SELECTED SO THAT THE COEFFICIENTS OF
C...  THE U1  TERMS SUM TO ONE (SINCE WE ARE INTERESTED IN OBTAINING
C...        X
C...  A FINITE DIFFERENCE APPROXIMATION FOR THIS FIRST DERIVATIVE).
C...  ALSO, WE SELECT A AND B SO THAT THE COEFFICIENTS OF THE U1
C...                                                            2X
C...  TERMS SUM TO ZERO IN ORDER TO DROP OUT THE CONTRIBUTION OF THIS
C...  SECOND DERIVATIVE (THE BASIC IDEA IS TO DROP OUT AS MANY OF THE
C...  DERIVATIVES AS POSSIBLE IN THE TAYLOR SERIES BEYOND THE DERI-
C...  VATIVE OF INTEREST, IN THIS CASE U1 , IN ORDER TO PRODUCE A
C...                                     X
C...  FINITE DIFFERENCE APPROXIMATION FOR THE DERIVATIVE OF MAXIMUM
C...  ACCURACY).  IN THIS CASE WE HAVE ONLY TWO CONSTANTS, A AND B,
C...  TO SELECT SO WE CAN DROP OUT ONLY THE SECOND DERIVATIVE, U1  ,
C...                                                             2X
C...  IN THE TAYLOR SERIES (IN ADDITION TO RETAINING THE FIRST DERI-
C...  VATIVE).  THIS PROCEDURE LEADS TO TWO LINEAR ALGEBRAIC EQUA-
C...  TIONS IN THE TWO CONSTANTS
C...
C...  A + 2B = 1
C...
C...  A + 4B = 0
C...
C...  SOLUTION OF THESE EQUATIONS FOR A AND B GIVES
C...
C...  A = 2, B = -1/2
C...
C...  SOLUTION OF EQUATIONS (9) AND (10) FOR U1  WITH THESE VALUES OF
C...  A AND B GIVES EQUATION (1)               X
C...
C...                                      2
C...  U1 = (1/2DX)(-3U1 + 4U2 - U3) + O(DX )                      (1)
C...    X
C...               2
C...  THE TERM O(DX ) INDICATES A PRINCIPAL ERROR TERM DUE TO TRUNCA-
C...                                                2
C...  TION OF THE TAYLOR SERIES WHICH IS OF ORDER DX .  THIS TERM IN
C...                    2
C...  FACT EQUALS U1  DX /3F, WHICH IS EASILY OBTAINED IN DERIVING
C...                3X
C...  EQUATION (1).
C...
C...  THIS SAME BASIC PROCEDURE CAN NOW BE APPLIED TO THE DERIVATION
C...  OF EQUATIONS (2) AND (3).
C...
C...     (2)  INTERIOR POINT I
C...
C...                                    2           3
C...  A(UI-1 = UI + UI (-DX) + UI  (-DX) + UI  (-DX) + ...)
C...                  X  1F      2X  2F      3X  3F
C...
C...                                    2           3
C...  B(UI+1 = UI + UI ( DX) + UI  ( DX) + UI  ( DX) + ...)
C...                  X  1F      2X  2F      3X  3F
C...
C...  -A + B = 1
C...
C...   A + B = 0
C...
C...  A = -1/2, B = 1/2
C...                                   2
C...  UI  = (1/2DX)(UI+1 - UI-1) + O(DX )                         (2)
C...    X
C...
C...     (3)  RIGHT END, POINT I = N
C...
C...                                      2            3
C...  A(UN-2 = UN + UN (-2DX) + UN  (-2DX) + UN  (-2DX) + ...)
C...                  X   1F      2X   2F      3X   3F
C...
C...                                      2            3
C...  B(UN-1 = UN + UN ( -DX) + UN  ( -DX) + UN  ( -DX) + ...)
C...                  X   1F      2X   2F      3X   3F
C...
C...  -2A - B = 1
C...
C...   4A + B = 0
C...
C...   A = -2, B = 1/2
C...                                          2
C...  UN  = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX )                  (3)
C...    X
C...
C...  THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1), (2) AND (3) CAN
C...  BE SUMMARIZED AS
C...
C...          -3   4  -1
C...
C...     1/2  -1   0   1
C...
C...           1  -4   3
C...
C...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 2, M =
C...  1, P = 0, 1, 2 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER-
C...  ENTIATION, MATH. GAZ., VOL. 25, 1941).
C...
C...  EQUATIONS (1), (2) AND (3) CAN NOW BE PROGRAMMED TO GENERATE
C...  THE DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF
C...                  X
C...  SUBROUTINE DSS002, RESPECTIVELY).
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R2FDX,      U,     UX,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT
      DX=(XU-XL)/DFLOAT(N-1)
      R2FDX=1.D+00/(2.D+00*DX)
      NM1=N-1
C...
C...  EQUATION (1) (NOTE - THE RHS OF THE FINITE DIFFERENCE APPROXI-
C...  TIONS, EQUATIONS (1), (2) AND (3) HAVE BEEN FORMATTED SO THAT
C...  THE NUMERICAL WEIGHTING COEFFICIENTS CAN BE MORE EASILY ASSOCI-
C...  ATED WITH THE BICKLEY MATRIX LISTED ABOVE)
      UX(1)=R2FDX*
     1(     -3.D+00*U(  1)     +4.D+00*U(  2)     -1.D+00*U(  3))
C...
C...  EQUATION (2)
      DO 1 I=2,NM1
      UX(I)=R2FDX*
     1(     -1.D+00*U(I-1)     +0.D+00*U(  I)     +1.D+00*U(I+1))
1     CONTINUE
C...
C...  EQUATION (3)
      UX(N)=R2FDX*
     1(      1.D+00*U(N-2)     -4.D+00*U(N-1)     +3.D+00*U(  N))
      RETURN
      END
*DECK DSS004
      SUBROUTINE DSS004(XL,XU,N,U,UX)
C...
C...  SUBROUTINE DSS004 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL
C...  FIVE-POINT, FOURTH-ORDER FINITE DIFFERENCE APPROXIMATIONS
C...
C...  ARGUMENT LIST
C...
C...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
C...
C...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
C...
C...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
C...             BOUNDARY POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
C...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
C...             TO BE COMPUTED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
C...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
C...             (OUTPUT)
C...
C...  THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR
C...  POLYNOMIALS) ARE GIVEN IN SUBROUTINE DSS002.
C...
C...  FIVE-POINT FORMULAS
C...
C...     (1)  LEFT END, POINT I = 1
C...
C...                                   2            3            4
C...  A(U2 = U1 + U1  ( DX) + U1  ( DX)  + U1  ( DX)  + U1  ( DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  ( DX)  + U1  ( DX)  + U1  ( DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  B(U3 = U1 + U1  (2DX) + U1  (2DX)  + U1  (2DX)  + U1  (2DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (2DX)  + U1  (2DX)  + U1  (2DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  C(U4 = U1 + U1  (3DX) + U1  (3DX)  + U1  (3DX)  + U1  (3DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (3DX)  + U1  (3DX)  + U1  (3DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  D(U5 = U1 + U1  (4DX) + U1  (4DX)  + U1  (4DX)  + U1  (4DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (4DX)  + U1  (4DX)  + U1  (4DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...  CONSTANTS A, B, C AND D ARE SELECTED SO THAT THE COEFFICIENTS
C...  OF THE U1  TERMS SUM TO ONE AND THE COEFFICIENTS OF THE U1  ,
C...           X                                                2X
C...  U1   AND U1   TERMS SUM TO ZERO
C...    3X       4X
C...
C...  A +   2B +   3C +   4D = 1
C...
C...  A +   4B +   9C +  16D = 0
C...
C...  A +   8B +  27C +  64D = 0
C...
C...  A +  16B +  81C + 256D = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
C...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
C...                                                            4X
C...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
C...               X
C...                                                         4
C...  U1  = (1/12DX)(-25U1 + 48U2 - 36U3 + 16U4 - 3U5) + O(DX )   (1)
C...    X
C...
C...     (2)  INTERIOR POINT, I = 2
C...
C...                                   2            3            4
C...  A(U1 = U2 + U2  (-DX) + U2  (-DX)  + U2  (-DX)  + U2  (-DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U2  (-DX)  + U2  (-DX)  + U2  (-DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  B(U3 = U2 + U2  ( DX) + U2  ( DX)  + U2  ( DX)  + U2  ( DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U2  ( DX)  + U2  ( DX)  + U2  ( DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  C(U4 = U2 + U2  (2DX) + U2  (2DX)  + U2  (2DX)  + U2  (2DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U2  (2DX)  + U2  (2DX)  + U2  (2DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                   2            3            4
C...  D(U5 = U2 + U2  (3DX) + U2  (3DX)  + U2  (3DX)  + U2  (3DX)
C...                X   1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U2  (3DX)  + U2  (3DX)  + U2  (3DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...  -A +   B +  2C +  3D = 1
C...
C...   A +   B +  4C +  9D = 0
C...
C...  -A +   B +  8C + 27D = 0
C...
C...   A +   B + 16C + 81D = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
C...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
C...                                                            4X
C...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
C...               X
C...                                                        4
C...  U2  = (1/12DX)(-3U1 - 10U2 + 18U3 -  6U4 +  U5) + O(DX )    (2)
C...    X
C...
C...     (3)  INTERIOR POINT I, I NE 2, N-1
C...
C...                                        2             3
C...  A(UI-2 = UI + UI  (-2DX)  + UI  (-2DX)  + UI  (-2DX)
C...                  X    1F       2X   2F       3X   3F
C...
C...                          4             5             6
C...              + UI  (-2DX)  + UI  (-2DX)  + UI  (-2DX)  + ...)
C...                  4X   4F       5X   5F       6X   6F
C...
C...                                        2             3
C...  B(UI-1 = UI + UI  ( -DX)  + UI  ( -DX)  + UI  ( -DX)
C...                  X    1F       2X   2F       3X   3F
C...
C...                          4             5             6
C...              + UI  ( -DX)  + UI  ( -DX)  + UI  ( -DX)  + ...)
C...                  4X   4F       5X   5F       6X   6F
C...
C...                                        2             3
C...  C(UI+1 = UI + UI  (  DX)  + UI  (  DX)  + UI  (  DX)
C...                  X    1F       2X   2F       3X   3F
C...
C...                          4             5             6
C...              + UI  (  DX)  + UI  (  DX)  + UI  (  DX)  + ...)
C...                  4X   4F       5X   5F       6X   6F
C...
C...                                        2             3
C...  D(UI+2 = UI + UI  ( 2DX)  + UI  ( 2DX)  + UI  ( 2DX)
C...                  X    1F       2X   2F       3X   3F
C...
C...                          4             5             6
C...              + UI  ( 2DX)  + UI  ( 2DX)  + UI  ( 2DX)  + ...)
C...                  4X   4F       5X   5F       6X   6F
C...
C...   -2A -   B +   C +  2D = 1
C...
C...    4A +   B +   C +  4D = 0
C...
C...   -8A -   B +   C +  8D = 0
C...
C...   16A +   B +   C + 16D = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
C...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
C...                                                            4X
C...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
C...               X
C...                                                          4
C...  UI  = (1/12DX)(UI-2 - 8UI-1 + 0UI + 8UI+1 - UI+2) + O(DX )  (3)
C...    X
C...
C...     (4)  INTERIOR POINT, I = N-1
C...
C...                                              2               3
C...  A(UN-4 = UN-1 + UN-1  (-3DX)  + UN-1  (-3DX)  + UN-1  (-3DX)
C...                      X    1F         2X   2F         3X   3F
C...
C...                       4               5               6
C...         + UN-1  (-3DX)  + UN-1  (-3DX)  + UN-1  (-3DX)  + ...
C...               4X   4F         5X   5F         6X   6F
C...
C...                                              2               3
C...  B(UN-3 = UN-1 + UN-1  (-2DX)  + UN-1  (-2DX)  + UN-1  (-2DX)
C...                      X    1F         2X   2F         3X   3F
C...
C...                       4               5               6
C...         + UN-1  (-2DX)  + UN-1  (-2DX)  + UN-1  (-2DX)  + ...
C...               4X   4F         5X   5F         6X   6F
C...
C...                                              2               3
C...  C(UN-2 = UN-1 + UN-1  ( -DX)  + UN-1  (- -X)  + UN-1  ( -DX)
C...                      X    1F         2X   2F         3X   3F
C...
C...                       4               5               6
C...         + UN-1  ( -DX)  + UN-1  ( -DX)  + UN-1  ( -DX)  + ...
C...               4X   4F         5X   5F         6X   6F
C...
C...                                              2               3
C...  D(UN   = UN-1 + UN-1  (  DX)  + UN-1  (  DX)  + UN-1  (  DX)
C...                      X    1F         2X   2F         3X   3F
C...
C...                       4               5               6
C...         + UN-1  (  DX)  + UN-1  (  DX)  + UN-1  (  DX)  + ...
C...               4X   4F         5X   5F         6X   6F
C...
C...  -3A -  2B -   C +   D = 1
C...
C...   9A +  4B +   C +   D = 0
C...
C... -27A -  8B -   C +   D = 0
C...
C...  81A + 16B +   C +   D = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
C...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
C...                                                            4X
C...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
C...               X
C...                                                                4
C...  UN-1  = (1/12DX)(-UN-4 + 6UN-3 - 18UN-2 + 10UN-1 + 3UN) + O(DX )
C...      X
C...                                                              (4)
C...
C...    (5)  RIGHT END, POINT I = N
C...
C...                                       2             3
C...  A(UN-4 = UN + UN (-4DX)  + UN  (-4DX)  + UN  (-4DX)
C...                  X   1F       2X   2F       3X   3F
C...
C...                         4             5             6
C...             + UN  (-4DX)  + UN  (-4DX)  + UN  (-4DX)  + ...)
C...                 4X   4F       5X   5F       6X   6F
C...
C...                                       2             3
C...  B(UN-3 = UN + UN (-3DX)  + UN  (-3DX)  + UN  (-3DX)
C...                  X   1F       2X   2F       3X   3F
C...
C...                         4             5             6
C...             + UN  (-3DX)  + UN  (-3DX)  + UN  (-3DX)  + ...)
C...                 4X   4F       5X   5F       6X   6F
C...
C...                                       2             3
C...  C(UN-2 = UN + UN (-2DX)  + UN  (-2DX)  + UN  (-2DX)
C...                  X   1F       2X   2F       3X   3F
C...
C...                         4             5             6
C...             + UN  (-2DX)  + UN  (-2DX)  + UN  (-2DX)  + ...)
C...                 4X   4F       5X   5F       6X   6F
C...
C...                                       2             3
C...  D(UN-1 = UN + UN ( -DX)  + UN  ( -DX)  + UN  ( -DX)
C...                  X   1F       2X   2F       3X   3F
C...
C...                         4             5             6
C...             + UN  ( -DX)  + UN  ( -DX)  + UN  ( -DX)  + ...)
C...                 4X   4F       5X   5F       6X   6F
C...
C...   -4A -  3B -  2C -   D = 1
C...
C...   16A +  9B +  4C +   D = 0
C...
C...  -64A - 27B -  8C -   D = 0
C...
C...  256A + 81B + 16C +   D = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU-
C...  TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U
C...                                                            4X
C...  TERMS, FOR U1  GIVES THE FOLLOWING FIVE-POINT APPROXIMATION
C...               X
C...                                                                4
C...  UN  = (1/12DX)(3UN-4 - 16UN-3 + 36UN-2 - 48UN-1 + 25UN) + O(DX )
C...    X
C...                                                              (5)
C...
C...  THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1) TO (5) CAN BE
C...  SUMMARIZED AS
C...
C...             -25   48  -36   16   -3
C...
C...              -3  -10   18   -6    1
C...
C...       1/12    1   -8    0    8   -1
C...
C...              -1    6  -18   10    3
C...
C...               3  -16   36  -48   25
C...
C...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 4, M =
C...  1, P = 0, 1, 2, 3, 4 (BICKLEY, W. G., FORMULAE FOR NUMERICAL
C...  DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941.  NOTE - THE BICKLEY
C...  COEFFICIENTS HAVE BEEN DIVIDED BY A COMMON FACTOR OF TWO).
C...
C...  EQUATIONS (1) TO (5) CAN NOW BE PROGRAMMED TO GENERATE THE
C...  DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF SUB-
C...              X
C...  ROUTINE DSS004 RESPECTIVELY).
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R4FDX,      U,     UX,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT
      DX=(XU-XL)/DFLOAT(N-1)
      R4FDX=1.D+00/(12.D+00*DX)
      NM2=N-2
C...
C...  EQUATION (1) (NOTE - THE RHS OF EQUATIONS (1), (2), (3), (4)
C...  AND (5) HAVE BEEN FORMATTED SO THAT THE NUMERICAL WEIGHTING
C...  COEFFICIENTS CAN BE MORE EASILY ASSOCIATED WITH THE BICKLEY
C...  MATRIX ABOVE)
      UX(  1)=R4FDX*
     1(-25.D+00      *U(  1)
     2 +48.D+00      *U(  2)
     3 -36.D+00      *U(  3)
     4 +16.D+00      *U(  4)
     5  -3.D+00      *U(  5))
C...
C...  EQUATION (2)
      UX(  2)=R4FDX*
     1( -3.D+00      *U(  1)
     2 -10.D+00      *U(  2)
     3 +18.D+00      *U(  3)
     4  -6.D+00      *U(  4)
     5  +1.D+00      *U(  5))
C...
C...  EQUATION (3)
      DO 1 I=3,NM2
      UX(  I)=R4FDX*
     1( +1.D+00      *U(I-2)
     2  -8.D+00      *U(I-1)
     3  +0.D+00      *U(  I)
     4  +8.D+00      *U(I+1)
     5  -1.D+00      *U(I+2))
1     CONTINUE
C...
C...  EQUATION (4)
      UX(N-1)=R4FDX*
     1( -1.D+00      *U(N-4)
     2  +6.D+00      *U(N-3)
     3 -18.D+00      *U(N-2)
     4 +10.D+00      *U(N-1)
     5  +3.D+00      *U(N  ))
C...
C...  EQUATION (5)
      UX(  N)=R4FDX*
     1( +3.D+00      *U(N-4)
     2 -16.D+00      *U(N-3)
     3 +36.D+00      *U(N-2)
     4 -48.D+00      *U(N-1)
     5 +25.D+00      *U(N  ))
      RETURN
      END
*DECK DSS006
      SUBROUTINE DSS006(XL,XU,N,U,UX)
C...
C...  SUBROUTINE DSS006 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL
C...  SEVEN-POINT, SIXTH-ORDER FINITE DIFFERENCE APPROXIMATIONS
C...
C...  ARGUMENT LIST
C...
C...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
C...
C...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
C...
C...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
C...             BOUNDARY POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
C...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
C...             TO BE COMPUTED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
C...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
C...             (OUTPUT)
C...
C...  THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR
C...  POLYNOMIALS) ARE GIVEN IN SUBROUTINES DSS002 AND DSS004.
C...
C...  SEVEN-POINT FORMULAS
C...
C...     (1)  LEFT END, POINT I = 1
C...
C...                                  2            3            4
C...  A(U2 = U1 + U1 ( DX) + U1  ( DX)  + U1  ( DX)  + U1  ( DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  ( DX)  + U1  ( DX)  + U1  ( DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  B(U3 = U1 + U1 (2DX) + U1  (2DX)  + U1  (2DX)  + U1  (2DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (2DX)  + U1  (2DX)  + U1  (2DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  C(U4 = U1 + U1 (3DX) + U1  (3DX)  + U1  (3DX)  + U1  (3DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (3DX)  + U1  (3DX)  + U1  (3DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  D(U5 = U1 + U1 (4DX) + U1  (4DX)  + U1  (4DX)  + U1  (4DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (4DX)  + U1  (4DX)  + U1  (4DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  E(U6 = U1 + U1 (5DX) + U1  (5DX)  + U1  (5DX)  + U1  (5DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (5DX)  + U1  (5DX)  + U1  (5DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  F(U7 = U1 + U1 (6DX) + U1  (6DX)  + U1  (6DX)  + U1  (6DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (6DX)  + U1  (6DX)  + U1  (6DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...  CONSTANTS A, B, C, D, E AND F ARE SELECTED SO THAT THE COEFFI-
C...  CIENTS  OF THE U1  TERMS SUM TO ONE AND THE COEFFICIENTS OF
C...                   X
C...  THE U1  , U1  , U1  , U1   AND U1   TERMS SUM TO ZERO.
C...        2X    3X    4X    5X       6X
C...
C...        1      1      1      1      1
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 1
C...
C...        2      2      2      2      2
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 0
C...
C...        3      3      3      3      3
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 0
C...
C...        4      4      4      4      4
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 0
C...
C...        5      5      5      5      5
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 0
C...
C...        6      6      6      6      6
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR U1  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X              X
C...  TION
C...
C...  U1  = (1/6F)(-1764*U1 + 4320U2 - 5400U3 + 4800U4
C...    X                                         6               (1)
C...              + 2700U5 + 864U6 - 120U7) + O(DX )
C...
C...     (2)  INTERIOR POINT, I = 2
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = 2
C...  AS
C...
C...  A(U1 = U2 + ...)
C...
C...  B(U3 = U2 + ...)
C...
C...  C(U4 = U2 + ...)
C...
C...  D(U5 = U2 + ...)
C...
C...  E(U6 = U2 + ...)
C...
C...  F(U7 = U2 + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 1
C...
C...    2      2      2      2      2      2
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 0
C...
C...    3      3      3      3      3      3
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 0
C...
C...    4      4      4      4      4      4
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 0
C...
C...    5      5      5      5      5      5
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 0
C...
C...    6      6      6      6      6      6
C...  -1 A +  1 B +  2 C +  3 D +  4 E +  5 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR U2  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X              X
C...  TION
C...
C...  U2  = (1/6F)(-120U1 - 924U2 + 1800U3 - 1200U4
C...    X                                       6                 (2)
C...              + 600U5 - 180U6 + 24U7) + O(DX )
C...
C...     (3)  INTERIOR POINT, I = 3
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = 3
C...  AS
C...
C...  A(U1 = U3 + ...)
C...
C...  B(U2 = U3 + ...)
C...
C...  C(U4 = U3 + ...)
C...
C...  D(U5 = U3 + ...)
C...
C...  E(U6 = U3 + ...)
C...
C...  F(U7 = U3 + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 1
C...
C...    2      2      2      2      2      2
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 0
C...
C...    3      3      3      3      3      3
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 0
C...
C...    4      4      4      4      4      4
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 0
C...
C...    5      5      5      5      5      5
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 0
C...
C...    6      6      6      6      6      6
C...  -2 A + -1 B +  1 C +  2 D +  3 E +  4 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR U3  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X              X
C...  TION
C...
C...  U3  = (1/6F)(24U1 - 288U2 - 420U3 + 960U4
C...    X                                     6                   (3)
C...             - 360U5 + 96U6 - 12U7) + O(DX )
C...
C...     (4)  INTERIOR POINT, I NE 2, 3, N-2, N-1
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = I
C...  AS
C...
C...  A(UI-3 = UI + ...)
C...
C...  B(UI-2 = UI + ...)
C...
C...  C(UI-1 = UI + ...)
C...
C...  D(UI+1 = UI + ...)
C...
C...  E(UI+2 = UI + ...)
C...
C...  F(UI+3 = UI + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 1
C...
C...    2      2      2      2      2      2
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 0
C...
C...    3      3      3      3      3      3
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 0
C...
C...    4      4      4      4      4      4
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 0
C...
C...    5      5      5      5      5      5
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 0
C...
C...    6      6      6      6      6      6
C...  -3 A + -2 B + -1 C +  1 D +  2 E +  3 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR UI  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X              X
C...  TION
C...
C...  UI  = (1/6F)(-12UI-3 + 108UI-2 - 540UI-1 + 0UI
C...    X                                             6           (4)
C...              + 540UI+1 - 108UI+2 + 12UI+3) + O(DX )
C...
C...     (5)  INTERIOR POINT, I = N-2
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N-2
C...  AS
C...
C...  A(UN-6 = UN-2 + ...)
C...
C...  B(UN-5 = UN-2 + ...)
C...
C...  C(UN-4 = UN-2 + ...)
C...
C...  D(UN-3 = UN-2 + ...)
C...
C...  E(UN-1 = UN-2 + ...)
C...
C...  F(UN   = UN-2 + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 1
C...
C...    2      2      2      2      2      2
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 0
C...
C...    3      3      3      3      3      3
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 0
C...
C...    4      4      4      4      4      4
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 0
C...
C...    5      5      5      5      5      5
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 0
C...
C...    6      6      6      6      6      6
C...  -4 A + -3 B + -2 C + -1 D +  1 E +  2 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR UN-2  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X                X
C...  TION
C...
C...  UN-2  = (1/6F)(12UN-6 - 96UN-5 + 360UN-4 - 960UN-3
C...      X                                          6            (5)
C...               + 420UN-2 + 288UN-1 - 24UN) + O(DX )
C...
C...     (6)  INTERIOR POINT, I = N-1
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N-1
C...  AS
C...
C...  A(UN-6 = UN-1 + ...)
C...
C...  B(UN-5 = UN-1 + ...)
C...
C...  C(UN-4 = UN-1 + ...)
C...
C...  D(UN-3 = UN-1 + ...)
C...
C...  E(UN-2 = UN-1 + ...)
C...
C...  F(UN   = UN-1 + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 1
C...
C...    2      2      2      2      2      2
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 0
C...
C...    3      3      3      3      3      3
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 0
C...
C...    4      4      4      4      4      4
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 0
C...
C...    5      5      5      5      5      5
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 0
C...
C...    6      6      6      6      6      6
C...  -5 A + -4 B + -3 C + -2 D + -1 E +  1 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR UN-1  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X                X
C...  TION
C...
C...  UN-1  = (1/6F)(-24UN-6 + 180UN-5 - 600UN-4 + 1200UN-3
C...      X                                             6         (6)
C...                - 1800UN-2 + 924UN-1 + 120UN) + O(DX )
C...
C...     (7)  RIGHT END, POINT I = N
C...
C...  THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N
C...  AS
C...
C...  A(UN-6 = UN   + ...)
C...
C...  B(UN-5 = UN   + ...)
C...
C...  C(UN-4 = UN   + ...)
C...
C...  D(UN-3 = UN   + ...)
C...
C...  E(UN-2 = UN   + ...)
C...
C...  F(UN-1 = UN   + ...)
C...
C...  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...    1      1      1      1      1      1
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 1
C...
C...    2      2      2      2      2      2
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0
C...
C...    3      3      3      3      3      3
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0
C...
C...    4      4      4      4      4      4
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0
C...
C...    5      5      5      5      5      5
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0
C...
C...    6      6      6      6      6      6
C...  -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE
C...  SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR UN  GIVES THE FOLLOWING SEVEN-POINT APPROXI-
C...   6X              X
C...  TION
C...
C...  UN  = (1/6F)(120UN-6 - 864UN-5 + 2700UN-4 - 4800UN-3
C...    X                                              6          (7)
C...             + 5400UN-2 - 4320UN-1 + 1764UN) + O(DX )
C...
C...  THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1) TO (7) CAN BE
C...  SUMMARIZED AS
C...
C...            -1764   4320  -5400   4800  -2700    864   -120
C...
C...             -120   -924   1800  -1200    600   -180     24
C...
C...               24   -288   -420    960   -360     96    -12
C...
C...     1/6F     -12    108   -540      0    540   -108     12
C...
C...               12    -96    360   -960    420    288    -24
C...
C...              -24    180   -600   1200  -1800    924    120
C...
C...              120   -864   2700  -4800   5400  -4320   1764
C...
C...  WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 6, M =
C...  1, P = 0 TO 6 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER-
C...  ENTIATION, MATH. GAZ., VOL. 25, 1941).
C...
C...  EQUATIONS (1) TO (7) CAN NOW BE PROGRAMMED TO GENERATE THE
C...  DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF SUB-
C...              X
C...  ROUTINE DSS006, RESPECTIVELY).
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R6FDX,      U,     UX,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT
      DX=(XU-XL)/DFLOAT(N-1)
      R6FDX=1.D+00/(720.D+00*DX)
      NM3=N-3
C...
C...  EQUATION (1)
      UX(  1)=R6FDX*
     1( -1764.D+00      *U(  1)
     2  +4320.D+00      *U(  2)
     3  -5400.D+00      *U(  3)
     4  +4800.D+00      *U(  4)
     5  -2700.D+00      *U(  5)
     6   +864.D+00      *U(  6)
     7   -120.D+00      *U(  7))
C...
C...  EQUATION (2)
      UX(  2)=R6FDX*
     1(  -120.D+00      *U(  1)
     2   -924.D+00      *U(  2)
     3  +1800.D+00      *U(  3)
     4  -1200.D+00      *U(  4)
     5   +600.D+00      *U(  5)
     6   -180.D+00      *U(  6)
     7    +24.D+00      *U(  7))
C...
C...  EQUATION (3)
      UX(  3)=R6FDX*
     1(   +24.D+00      *U(  1)
     2   -288.D+00      *U(  2)
     3   -420.D+00      *U(  3)
     4   +960.D+00      *U(  4)
     5   -360.D+00      *U(  5)
     6    +96.D+00      *U(  6)
     7    -12.D+00      *U(  7))
C...
C...  EQUATION (4)
      DO 1 I=4,NM3
      UX(  I)=R6FDX*
     1(   -12.D+00      *U(I-3)
     2   +108.D+00      *U(I-2)
     3   -540.D+00      *U(I-1)
     4     +0.D+00      *U(I  )
     5   +540.D+00      *U(I+1)
     6   -108.D+00      *U(I+2)
     7    +12.D+00      *U(I+3))
1     CONTINUE
C...
C...  EQUATION (5)
      UX(N-2)=R6FDX*
     1(   +12.D+00      *U(N-6)
     2    -96.D+00      *U(N-5)
     3   +360.D+00      *U(N-4)
     4   -960.D+00      *U(N-3)
     5   +420.D+00      *U(N-2)
     6   +288.D+00      *U(N-1)
     7    -24.D+00      *U(N  ))
C...
C...  EQUATION (6)
      UX(N-1)=R6FDX*
     1(   -24.D+00      *U(N-6)
     2   +180.D+00      *U(N-5)
     3   -600.D+00      *U(N-4)
     4  +1200.D+00      *U(N-3)
     5  -1800.D+00      *U(N-2)
     6   +924.D+00      *U(N-1)
     7   +120.D+00      *U(N  ))
C...
C...  EQUATION (7)
      UX(  N)=R6FDX*
     1(  +120.D+00      *U(N-6)
     2   -864.D+00      *U(N-5)
     3  +2700.D+00      *U(N-4)
     4  -4800.D+00      *U(N-3)
     5  +5400.D+00      *U(N-2)
     6  -4320.D+00      *U(N-1)
     7  +1764.D+00      *U(N  ))
      RETURN
      END
*DECK DSS008
      SUBROUTINE DSS008(XL,XU,N,U,UX)
C...
C...  SUBROUTINE DSS008 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL
C...  NINE-POINT, EIGHTH-ORDER FINITE DIFFERENCE APPROXIMATIONS
C...
C...  ARGUMENT LIST
C...
C...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
C...
C...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
C...
C...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
C...             BOUNDARY POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
C...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
C...             TO BE COMPUTED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
C...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
C...             (OUTPUT)
C...
C...  THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR
C...  POLYNOMIALS) ARE GIVEN IN SUBROUTINES DSS002, DSS004 AND
C...  DSS006.
C...
C...  NINE-POINT FORMULAS
C...
C...     (1)  LEFT END, POINT I = 1
C...
C...                                  2            3            4
C...  A(U2 = U1 + U1 ( DX) + U1  ( DX)  + U1  ( DX)  + U1  ( DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  ( DX)  + U1  ( DX)  + U1  ( DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  B(U3 = U1 + U1 (2DX) + U1  (2DX)  + U1  (2DX)  + U1  (2DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (2DX)  + U1  (2DX)  + U1  (2DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  C(U4 = U1 + U1 (3DX) + U1  (3DX)  + U1  (3DX)  + U1  (3DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (3DX)  + U1  (3DX)  + U1  (3DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  D(U5 = U1 + U1 (4DX) + U1  (4DX)  + U1  (4DX)  + U1  (4DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (4DX)  + U1  (4DX)  + U1  (4DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  E(U6 = U1 + U1 (5DX) + U1  (5DX)  + U1  (5DX)  + U1  (5DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (5DX)  + U1  (5DX)  + U1  (5DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...                                  2            3            4
C...  F(U7 = U1 + U1 (6DX) + U1  (6DX)  + U1  (6DX)  + U1  (6DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (6DX)  + U1  (6DX)  + U1  (6DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  G(U8 = U1 + U1 (7DX) + U1  (7DX)  + U1  (7DX)  + U1  (7DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (7DX)  + U1  (7DX)  + U1  (7DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...                                  2            3            4
C...  H(U9 = U1 + U1 (8DX) + U1  (8DX)  + U1  (8DX)  + U1  (8DX)
C...                X  1F      2X  2F       3X  3F       4X  4F
C...
C...                      5            6            7
C...           + U1  (8DX)  + U1  (8DX)  + U1  (8DX)  + ...)
C...               5X  5F       6X  6F       7X  7F
C...
C...  CONSTANTS A, B, C, D, E, F, G AND H ARE SELECTED SO THAT THE
C...  COEFFICIENTS OF THE U1  TERMS SUM TO ONE AND THE COEFFICIENTS OF
C...                        X
C...  THE U1  , U1  , U1  , U1  , U1  , U1   AND U1   SUM TO ZERO.
C...        2X    3X    4X    5X    6X,   7X       8X
C...
C...        1      1      1      1      1      1      1
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 1
C...
C...        2      2      2      2      2      2      2
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        3      3      3      3      3      3      3
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        4      4      4      4      4      4      4
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        5      5      5      5      5      5      5
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        6      6      6      6      6      6      6
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        7      7      7      7      7      7      7
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...        8      8      8      8      8      8      8
C...  A +  2 B +  3 C +  4 D +  5 E +  6 F +  7 G +  8 H = 0
C...
C...  SIMULTANEOUS SOLUTION FOR A, B, C, D, E, F, G AND H FOLLOWED BY
C...  THE SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE
C...  U   TERMS, FOR U1  GIVES THE FOLLOWING NINE-POINT APPROXIMATION
C...   8X              X
C...
C...  U1  = (1/8F)(-109584U1 + 322560U2 - 564480U3 + 752640U4
C...    X
C...               -705600U5 + 451584U6 - 188160U7 + 46080U8
C...                               8
C...               - 5040U9) + O(DX )                             (1)
C...
C...  THE PRECEDING ANALYSIS CAN BE REPEATED TO PRODUCE NINE-POINT
C...  APPROXIMATIONS FOR THE FIRST DERIVATIVES U2 , U3 , U4 , UI ,
C...                                             X    X    X    X
C...  UN-3 , UN-2 , UN-1  AND UN .  THE RESULTS CAN BE SUMMARIZED BY
C...      X      X      X       X
C...  THE FOLLOWING BICKLEY MATRIX FOR N = 8, M = 1, P = 0 TO 8,
C...  (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFERENTIATION, MATH.
C...  GAZ., VOL. 25, 1941)
C...
C...            -109584   322560  -564480   752640  -705600
C...
C...              -5040   -64224   141120  -141120   117600
C...
C...                720   -11520   -38304    80640   -50400
C...
C...               -240     2880   -20160   -18144    50400
C...
C...     1/8F       144    -1536     8064   -32256        0
C...
C...               -144     1440    -6720    20160   -50400
C...
C...                240    -2304    10080   -26880    50400
C...
C...               -720     6720   -28224    70560  -117600
C...
C...               5040   -46080   188160  -451584   705600
C...
C...  FROM THIS BICKLEY MATRIX, THE FINITE DIFFERENCE APPROXIMATION OF
C...  THE FIRST DERIVATIVE CAN BE PROGRAMMED FOR EACH OF THE GRID POINTS
C...  1, 2, 3, 4,..., I,..., N-3, N-2, N-1, N (TAKING INTO ACCOUNT
C...  THE SYMMETRY PROPERTIES OF THE BICKLEY MATRIX).
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R8FDX,      U,     UX,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT
      DX=(XU-XL)/DFLOAT(N-1)
      R8FDX=1.D+00/(40320.D+00*DX)
      NM4=N-4
C...
C...  GRID POINT 1
      UX(  1)=R8FDX*
     1   (-109584.D+00     *U(  1)
     2    +322560.D+00     *U(  2)
     3    -564480.D+00     *U(  3)
     4    +752640.D+00     *U(  4)
     5    -705600.D+00     *U(  5)
     6    +451584.D+00     *U(  6)
     7    -188160.D+00     *U(  7)
     8    +46080.D+00      *U(  8)
     9    -5040.D+00       *U(  9))
C...
C...  GRID POINT 2
      UX(  2)=R8FDX*
     1   (-5040.D+00       *U(  1)
     2    -64224.D+00      *U(  2)
     3    +141120.D+00     *U(  3)
     4    -141120.D+00     *U(  4)
     5    +117600.D+00     *U(  5)
     6    -70560.D+00      *U(  6)
     7    +28224.D+00      *U(  7)
     8    -6720.D+00       *U(  8)
     9    +720.D+00        *U(  9))
C...
C...  GRID POINT 3
      UX(  3)=R8FDX*
     1   (+720.D+00        *U(  1)
     2    -11520.D+00      *U(  2)
     3    -38304.D+00      *U(  3)
     4    +80640.D+00      *U(  4)
     5    -50400.D+00      *U(  5)
     6    +26880.D+00      *U(  6)
     7    -10080.D+00      *U(  7)
     8    +2304.D+00       *U(  8)
     9    -240.D+00        *U(  9))
C...
C...  GRID POINT 4
      UX(  4)=R8FDX*
     1   (-240.D+00        *U(  1)
     2    +2880.D+00       *U(  2)
     3    -20160.D+00      *U(  3)
     4    -18144.D+00      *U(  4)
     5    +50400.D+00      *U(  5)
     6    -20160.D+00      *U(  6)
     7    +6720.D+00       *U(  7)
     8    -1440.D+00       *U(  8)
     9    +144.D+00        *U(  9))
C...
C...  GRID POINT I
      DO 1 I=5,NM4
      UX(  I)=R8FDX*
     1   (+144.D+00        *U(I-4)
     2    -1536.D+00       *U(I-3)
     3    +8064.D+00       *U(I-2)
     4    -32256.D+00      *U(I-1)
     5    +0.D+00          *U(I  )
     6    +32256.D+00      *U(I+1)
     7    -8064.D+00       *U(I+2)
     8    +1536.D+00       *U(I+3)
     9    -144.D+00        *U(I+4))
1     CONTINUE
C...
C...  GRID POINT N-3
      UX(N-3)=R8FDX*
     1   (-144.D+00        *U(N-8)
     2    +1440.D+00       *U(N-7)
     3    -6720.D+00       *U(N-6)
     4    +20160.D+00      *U(N-5)
     5    -50400.D+00      *U(N-4)
     6    +18144.D+00      *U(N-3)
     7    +20160.D+00      *U(N-2)
     8    -2880.D+00       *U(N-1)
     9    +240.D+00        *U(N  ))
C...
C...  GRID POINT N-2
      UX(N-2)=R8FDX*
     1   (+240.D+00        *U(N-8)
     2    -2304.D+00       *U(N-7)
     3    +10080.D+00      *U(N-6)
     4    -26880.D+00      *U(N-5)
     5    +50400.D+00      *U(N-4)
     6    -80640.D+00      *U(N-3)
     7    +38304.D+00      *U(N-2)
     8    +11520.D+00      *U(N-1)
     9    -720.D+00        *U(N  ))
C...
C...  GRID POINT N-1
      UX(N-1)=R8FDX*
     1   (-720.D+00        *U(N-8)
     2    +6720.D+00       *U(N-7)
     3    -28224.D+00      *U(N-6)
     4    +70560.D+00      *U(N-5)
     5    -117600.D+00     *U(N-4)
     6    +141120.D+00     *U(N-3)
     7    -141120.D+00     *U(N-2)
     8    +64224.D+00      *U(N-1)
     9    +5040.D+00       *U(N  ))
C...
C...  GRID POINT N
      UX(N  )=R8FDX*
     1   (+5040.D+00       *U(N-8)
     2    -46080.D+00      *U(N-7)
     3    +188160.D+00     *U(N-6)
     4    -451584.D+00     *U(N-5)
     5    +705600.D+00     *U(N-4)
     6    -752640.D+00     *U(N-3)
     7    +564480.D+00     *U(N-2)
     8    -322560.D+00     *U(N-1)
     9    +109584.D+00     *U(N  ))
      RETURN
      END
*DECK DSS010
      SUBROUTINE DSS010(XL,XU,N,U,UX)
C...
C...  SUBROUTINE DSS010 COMPUTES THE FIRST DERIVATIVE, U , OF A
C...                                                    X
C...  VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL
C...  ELEVEN-POINT, TENTH-ORDER FINITE DIFFERENCE APPROXIMATIONS
C...
C...  ARGUMENT LIST
C...
C...     XL      LOWER BOUNDARY VALUE OF X (INPUT)
C...
C...     XU      UPPER BOUNDARY VALUE OF X (INPUT)
C...
C...     N       NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE
C...             BOUNDARY POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT
C...             THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS
C...             TO BE COMPUTED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL
C...             VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS
C...             (OUTPUT)
C...
C...  THE MATHEMATICAL DETAILS OF THE FINITE DIFFERENCE APPROXIMA-
C...  TIONS ARE GIVEN IN SUBROUTINES DSS002, DSS004, DSS006 AND
C...  DSS008, AND CAN BE SUMMARIZED BY THE FOLLOWING BICKLEY MATRIX
C...  FOR N = 10, M = 1, P = 0 TO 10 (BICKLEY, W. G., FORMULAE FOR
C...  NUMERICAL DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941)
C...
C...      -10626840            36288000           -81648000
C...      145152000          -190512000           182891520
C...
C...        -362880            -6636960            16329600
C...      -21772800            25401600           -22861440
C...
C...          40320             -806400            -4419360
C...        9676800            -8467200             6773760
C...
C...         -10080              151200            -1360800
C...       -2756160             6350400            -3810240
C...
C...           4320              -57600              388800
C...       -2073600            -1330560             4354560
C...
C...          -2880               36000             -216000
C...         864000            -3024000                   0
C...
C...           2880              -34560              194400
C...        -691200             1814400            -4354560
C...
C...          -4320               50400             -272160
C...         907200            -2116800             3810240
C...
C...          10080             -115200              604800
C...       -1935360             4233600            -6773760
C...
C...         -40320              453600            -2332800
C...        7257600           -15240960            22861440
C...
C...         362880            -4032000            20412000
C...      -62208000           127008000          -182891520
C...
C...  EACH ENTRY IN THIS TABLE SHOULD BE MULTIPLIED BY 1/10F TO
C...  OBTAIN THE FINAL WEIGHTING COEFFICIENTS.  FROM THIS BICKLEY
C...  MATRIX, THE FINITE DIFFERENCE APPROXIMATION OF THE FIRST
C...  DERIVATIVE CAN BE PROGRAMMED FOR EACH OF THE GRID POINTS 1, 2,
C...  3, 4, 5,..., I,..., N-4, N-3, N-2, N-1 AND N (TAKING INTO
C...  ACCOUNT THE SYMMETRY PROPERTIES OF THE MATRIX).
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX, R10FDX,      U,     UX,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT
      DX=(XU-XL)/DFLOAT(N-1)
      R10FDX=1.D+00/(3628800.D+00*DX)
      NM5=N-5
C...
C...  GRID POINT 1
      UX(  1)=R10FDX*
     1   (-10628640.D+00      *U(   1)
     2    +36288000.D+00      *U(   2)
     3    -81648000.D+00      *U(   3)
     4    +145152000.D+00     *U(   4)
     5    -190512000.D+00     *U(   5)
     6    +182891520.D+00     *U(   6)
     7    -127008000.D+00     *U(   7)
     8    +62208000.D+00      *U(   8)
     9    -20412000.D+00      *U(   9)
     A    +4032000.D+00       *U(  10)
     B    -362880.D+00        *U(  11))
C...
C...  GRID POINT 2
      UX(  2)=R10FDX*
     1   (-362880.D+00        *U(   1)
     2    -6636960.D+00       *U(   2)
     3    +16329600.D+00      *U(   3)
     4    -21772800.D+00      *U(   4)
     5    +25401600.D+00      *U(   5)
     6    -22861440.D+00      *U(   6)
     7    +15240960.D+00      *U(   7)
     8    -7257600.D+00       *U(   8)
     9    +2332800.D+00       *U(   9)
     A    -453600.D+00        *U(  10)
     B    +40320.D+00         *U(  11))
C...
C...  GRID POINT 3
      UX(  3)=R10FDX*
     1   (+40320.D+00         *U(   1)
     2    -806400.D+00        *U(   2)
     3    -4419360.D+00       *U(   3)
     4    +9676800.D+00       *U(   4)
     5    -8467200.D+00       *U(   5)
     6    +6773760.D+00       *U(   6)
     7    -4233600.D+00       *U(   7)
     8    +1935360.D+00       *U(   8)
     9    -604800.D+00        *U(   9)
     A    +115200.D+00        *U(  10)
     B    -10080.D+00         *U(  11))
C...
C...  GRID POINT 4
      UX(  4)=R10FDX*
     1   (-10080.D+00         *U(   1)
     2    +151200.D+00        *U(   2)
     3    -1360800.D+00       *U(   3)
     4    -2756160.D+00       *U(   4)
     5    +6350400.D+00       *U(   5)
     6    -3810240.D+00       *U(   6)
     7    +2116800.D+00       *U(   7)
     8    -907200.D+00        *U(   8)
     9    +272160.D+00        *U(   9)
     A    -50400.D+00         *U(  10)
     B    +4320.D+00          *U(  11))
C...
C...  GRID POINT 5
      UX(  5)=R10FDX*
     1   (+4320.D+00          *U(   1)
     2    -57600.D+00         *U(   2)
     3    +388800.D+00        *U(   3)
     4    -2073600.D+00       *U(   4)
     5    -1330560.D+00       *U(   5)
     6    +4354560.D+00       *U(   6)
     7    -1814400.D+00       *U(   7)
     8    +691200.D+00        *U(   8)
     9    -194400.D+00        *U(   9)
     A    +34560.D+00         *U(  10)
     B    -2880.D+00          *U(  11))
C...
C...  GRID POINT I, I NE 1 TO 5, N-4 TO N
      DO 1 I=6,NM5
      UX(  I)=R10FDX*
     1   (-2880.D+00          *U( I-5)
     2    +36000.D+00         *U( I-4)
     3    -216000.D+00        *U( I-3)
     4    +864000.D+00        *U( I-2)
     5    -3024000.D+00       *U( I-1)
     6    +0.D+00             *U( I  )
     7    +3024000.D+00       *U( I+1)
     8    -864000.D+00        *U( I+2)
     9    +216000.D+00        *U( I+3)
     A    -36000.D+00         *U( I+4)
     B    +2880.D+00          *U( I+5))
1     CONTINUE
C...
C...  GRID POINT N-4
      UX(N-4)=R10FDX*
     1   (+2880.D+00          *U(N-10)
     2    -34560.D+00         *U( N-9)
     3    +194400.D+00        *U( N-8)
     4    -691200.D+00        *U( N-7)
     5    +1814400.D+00       *U( N-6)
     6    -4354560.D+00       *U( N-5)
     7    +1330560.D+00       *U( N-4)
     8    +2073600.D+00       *U( N-3)
     9    -388800.D+00        *U( N-2)
     A    +57600.D+00         *U( N-1)
     B    -4320.D+00          *U( N  ))
C...
C...  GRID POINT N-3
      UX(N-3)=R10FDX*
     1   (-4320.D+00          *U(N-10)
     2    +50400.D+00         *U( N-9)
     3    -272160.D+00        *U( N-8)
     4    +907200.D+00        *U( N-7)
     5    -2116800.D+00       *U( N-6)
     6    +3810240.D+00       *U( N-5)
     7    -6350400.D+00       *U( N-4)
     8    +2756160.D+00       *U( N-3)
     9    +1360800.D+00       *U( N-2)
     A    -151200.D+00        *U( N-1)
     B    +10080.D+00         *U( N  ))
C...
C...  GRID POINT N-2
      UX(N-2)=R10FDX*
     1   (+10080.D+00         *U(N-10)
     2    -115200.D+00        *U( N-9)
     3    +604800.D+00        *U( N-8)
     4    -1935360.D+00       *U( N-7)
     5    +4233600.D+00       *U( N-6)
     6    -6773760.D+00       *U( N-5)
     7    +8467200.D+00       *U( N-4)
     8    -9676800.D+00       *U( N-3)
     9    +4419360.D+00       *U( N-2)
     A    +806400.D+00        *U( N-1)
     B    -40320.D+00         *U( N  ))
C...
C...  GRID POINT N-1
      UX(N-1)=R10FDX*
     1   (-40320.D+00         *U(N-10)
     2    +453600.D+00        *U( N-9)
     3    -2332800.D+00       *U( N-8)
     4    +7257600.D+00       *U( N-7)
     5    -15240960.D+00      *U( N-6)
     6    +22861440.D+00      *U( N-5)
     7    -25401600.D+00      *U( N-4)
     8    +21772800.D+00      *U( N-3)
     9    -16329600.D+00      *U( N-2)
     A    +6636960.D+00       *U( N-1)
     B    +362880.D+00        *U( N  ))
C...
C...  GRID POINT N
      UX(N  )=R10FDX*
     1   (+362880.D+00        *U(N-10)
     2    -4032000.D+00       *U( N-9)
     3    +20412000.D+00      *U( N-8)
     4    -62208000.D+00      *U( N-7)
     5    +127008000.D+00     *U( N-6)
     6    -182891520.D+00     *U( N-5)
     7    +190512000.D+00     *U( N-4)
     8    -145152000.D+00     *U( N-3)
     9    +81648000.D+00      *U( N-2)
     A    -36288000.D+00      *U( N-1)
     B    +10628640.D+00      *U( N  ))
      RETURN
      END
*DECK DSS012
      SUBROUTINE DSS012(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS012 IS AN APPLICATION OF FIRST-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS WITH THE
C...  SIMPLEST FORM
C...
C...                            U  + V*U  = 0                        (1)
C...                             T      X
C...
C...  THE FIRST FIVE PARAMETERS, XL, XU, N, U AND UX, ARE THE SAME
C...  AS FOR SUBROUTINES DSS002 TO DSS010 AS DEFINED IN THOSE ROUTINES.
C...  THE SIXTH PARAMETER, V, MUST BE PROVIDED TO DSS012 SO THAT THE
C...  DIRECTION OF FLOW IN EQUATION (1) CAN BE USED TO SELECT THE
C...  APPROPRIATE FINITE DIFFERENCE APPROXIMATION FOR THE FIRST-ORDER
C...  SPATIAL DERIVATIVE IN EQUATION (1), U .  THE CONVENTION FOR THE
C...  SIGN OF V IS                         X
C...
C...     FLOW LEFT TO RIGHT                 V GT 0
C...     (I.E., IN THE DIRECTION            (I.E., THE SIXTH ARGUMENT IS
C...     OF INCREASING X)                   POSITIVE IN CALLING DSS012)
C...
C...     FLOW RIGHT TO LEFT                 V LT 0
C...     (I.E., IN THE DIRECTION            (I.E., THE SIXTH ARGUMENT IS
C...     OF DECREASING X)                   NEGATIVE IN CALLING DSS012)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,      U,     UX,      V,     XL,     XU,
     1                 DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFERENCE
C...  APPROXIMATION DEPENDING ON THE SIGN OF V IN EQUATION (1).  THE
C...  ORIGIN OF THE FINITE DIFFERENCE APPROXIMATIONS USED BELOW IS GIVEN
C...  AT THE END OF SUBROUTINE DSS012.
      DX=(XU-XL)/DFLOAT(N-1)
      IF(V.LT.0.D+00)GO TO 10
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
              UX(1)=(U(2)-U(1))/DX
              DO 1 I=2,N
              UX(I)=(U(I)-U(I-1))/DX
1             CONTINUE
              RETURN
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
10            NM1=N-1
              DO 2 I=1,NM1
              UX(I)=(U(I+1)-U(I))/DX
2             CONTINUE
              UX(N)=(U(N)-U(N-1))/DX
              RETURN
C...
C...  THE BACKWARD DIFFERENCES IN SECTION (1) ABOVE ARE BASED ON THE
C...  TAYLOR SERIES
C...
C...                                  2           3
C...  UI-1 = UI + UI (-DX) + UI  (-DX) + UI  (-DX) + ...
C...                X  1F      2X  2F      3X  3F
C...
C...                                          2
C...  IF THIS SERIES IS TRUNCATED AFTER THE DX  TERM AND THE RESULTING
C...  EQUATION SOLVED FOR U ,  WE OBTAIN IMMEDIATELY
C...                       X
C...
C...  UI  = (UI - UI-1)/DX + O(DX)
C...    X
C...
C...  WHICH IS THE FIRST-ORDER BACKWARD DIFFERENCE USED IN DO LOOP 1.
C...  THE DERIVATIVE U1  IS COMPUTED BY USING THE POINT TO THE RIGHT OF
C...                   X
C...  U1, I.E., U2, SINCE THIS IS THE ONLY POINT AVAILABLE IF FICTITIOUS
C...  POINTS TO THE LEFT OF U1 ARE TO BE AVOIDED.
C...
C...  THE FORWARD DIFFERENCES IN SECTION (2) ABOVE ARE BASED ON THE
C...  TAYLOR SERIES
C...
C...                                  2           3
C...  UI+1 = UI + UI ( DX) + UI  ( DX) + UI  ( DX) + ...
C...                X  1F      2X  2F      3X  3F
C...
C...                                          2
C...  IF THIS SERIES IS TRUNCATED AFTER THE DX  TERM AND THE RESULTING
C...  EQUATION SOLVED FOR U ,  WE OBTAIN IMMEDIATELY
C...                       X
C...
C...  UI  = (UI+1 - UI)/DX + O(DX)
C...    X
C...
C...  WHICH IS THE FIRST-ORDER FORWARD DIFFERENCE USED IN DO LOOP 2.
C...  THE DERIVATIVE UN  IS COMPUTED BY USING THE POINT TO THE LEFT OF
C...                   X
C...  UN (UN-1), SINCE THIS IS THE ONLY POINT AVAILABLE IF FICTITIOUS
C...  POINTS TO THE RIGHT OF UN ARE TO BE AVOIDED.
      END
*DECK DSS014
      SUBROUTINE DSS014(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS014 IS AN APPLICATION OF SECOND-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
C...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
C...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
C...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
C...  GAZETTE, PP. 19-27, 1941, N = 2, M = 1, P = 0, 1, 2.
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R2FDX,      U,     UX,      V,     XL,
     1                     XU, DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
C...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
C...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      DX=(XU-XL)/DFLOAT(N-1)
      R2FDX=1.D+00/(2.D+00*DX)
      IF(V.LT.0.D+00)GO TO 10
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
      UX(1)=R2FDX*
     1(     -3.D+00     *U(  1)
     2      +4.D+00     *U(  2)
     3      -1.D+00     *U(  3))
      UX(2)=R2FDX*
     1(     -1.D+00     *U(  1)
     2      +0.D+00     *U(  2)
     3      +1.D+00     *U(  3))
      DO 1 I=3,N
      UX(I)=R2FDX*
     1(     +1.D+00     *U(I-2)
     2      -4.D+00     *U(I-1)
     3      +3.D+00     *U(I  ))
1     CONTINUE
      RETURN
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
10    NM2=N-2
      DO 2 I=1,NM2
      UX(I)=R2FDX*
     1(     -3.D+00     *U(I  )
     2      +4.D+00     *U(I+1)
     3      -1.D+00     *U(I+2))
2     CONTINUE
      UX(N-1)=R2FDX*
     1(     -1.D+00     *U(N-2)
     2      +0.D+00     *U(N-1)
     3      +1.D+00     *U(N  ))
      UX(N)=R2FDX*
     1(     +1.D+00     *U(N-2)
     2      -4.D+00     *U(N-1)
     3      +3.D+00     *U(N  ))
      RETURN
      END
*DECK DSS016
      SUBROUTINE DSS016(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS016 IS AN APPLICATION OF THIRD-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
C...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
C...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
C...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
C...  GAZETTE, PP. 19-27, 1941, N = 3, M = 1, P = 0, 1, 2, 3.
C...
      DOUBLE PRECISION     DX,  R3FDX,      U,     UX,      V,     XL,
     1                     XU, DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
C...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
C...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      DX=(XU-XL)/DFLOAT(N-1)
      R3FDX=1.D+00/(6.D+00*DX)
      IF(V.LT.0.D+00)GO TO 10
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
      UX(  1)=R3FDX*
     1( -11.D+00     *U(  1)
     2  +18.D+00     *U(  2)
     3   -9.D+00     *U(  3)
     4   +2.D+00     *U(  4))
      UX(  2)=R3FDX*
     1(  -2.D+00     *U(  1)
     2   -3.D+00     *U(  2)
     3   +6.D+00     *U(  3)
     4   -1.D+00     *U(  4))
      UX(  3)=R3FDX*
     1(  +1.D+00     *U(  1)
     2   -6.D+00     *U(  2)
     3   +3.D+00     *U(  3)
     4   +2.D+00     *U(  4))
      DO 1 I=4,N
      UX(  I)=R3FDX*
     1(  -2.D+00     *U(I-3)
     2   +9.D+00     *U(I-2)
     3  -18.D+00     *U(I-1)
     4  +11.D+00     *U(I  ))
1     CONTINUE
      RETURN
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
10    NM3=N-3
      DO 2 I=1,NM3
      UX(  I)=R3FDX*
     1( -11.D+00     *U(I  )
     2  +18.D+00     *U(I+1)
     3   -9.D+00     *U(I+2)
     4   +2.D+00     *U(I+3))
2     CONTINUE
      UX(N-2)=R3FDX*
     1(  -2.D+00     *U(N-3)
     2   -3.D+00     *U(N-2)
     3   +6.D+00     *U(N-1)
     4   -1.D+00     *U(N  ))
      UX(N-1)=R3FDX*
     1(  +1.D+00     *U(N-3)
     2   -6.D+00     *U(N-2)
     3   +3.D+00     *U(N-1)
     4   +2.D+00     *U(N  ))
      UX(  N)=R3FDX*
     1(  -2.D+00     *U(N-3)
     2   +9.D+00     *U(N-2)
     3  -18.D+00     *U(N-1)
     4  +11.D+00     *U(N  ))
      RETURN
      END
*DECK DSS018
      SUBROUTINE DSS018(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS018 IS AN APPLICATION OF THIRD-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
C...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
C...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
C...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
C...  GAZETTE, PP. 19-27, 1941, N = 3, M = 1, P = 0, 1, 2, 3.  THE
C...  IMPLEMENTATION IS THE **FOUR-POINT BIASED UPWIND FORMULA** OF
C...  M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE
C...  ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69,
C...  AUGUST, 1978
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION     DX,  R3FDX,      U,     UX,      V,     XL,
     1                     XU, DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
C...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
C...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      DX=(XU-XL)/DFLOAT(N-1)
      R3FDX=1.D+00/(6.D+00*DX)
      IF(V.LT.0.D+00)GO TO 10
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
      UX(  1)=R3FDX*
     1( -11.D+00     *U(  1)
     2  +18.D+00     *U(  2)
     3   -9.D+00     *U(  3)
     4   +2.D+00     *U(  4))
      UX(  2)=R3FDX*
     1(  -2.D+00     *U(  1)
     2   -3.D+00     *U(  2)
     3   +6.D+00     *U(  3)
     4   -1.D+00     *U(  4))
      NM1=N-1
      DO 1 I=3,NM1
      UX(  I)=R3FDX*
     1(  +1.D+00     *U(I-2)
     2   -6.D+00     *U(I-1)
     3   +3.D+00     *U(I  )
     4   +2.D+00     *U(I+1))
1     CONTINUE
      UX(  N)=R3FDX*
     1(  -2.D+00     *U(N-3)
     2   +9.D+00     *U(N-2)
     3  -18.D+00     *U(N-1)
     4  +11.D+00     *U(N  ))
      RETURN
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
10    UX(  1)=R3FDX*
     1( -11.D+00     *U(  1)
     2  +18.D+00     *U(  2)
     3   -9.D+00     *U(  3)
     4   +2.D+00     *U(  4))
      NM2=N-2
      DO 2 I=2,NM2
      UX(  I)=R3FDX*
     1(  -2.D+00     *U(I-1)
     2   -3.D+00     *U(I  )
     3   +6.D+00     *U(I+1)
     4   -1.D+00     *U(I+2))
2     CONTINUE
      UX(N-1)=R3FDX*
     1(  +1.D+00     *U(N-3)
     2   -6.D+00     *U(N-2)
     3   +3.D+00     *U(N-1)
     4   +2.D+00     *U(N  ))
      UX(  N)=R3FDX*
     1(  -2.D+00     *U(N-3)
     2   +9.D+00     *U(N-2)
     3  -18.D+00     *U(N-1)
     4  +11.D+00     *U(N  ))
      RETURN
      END
*DECK DSS020
      SUBROUTINE DSS020(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS020 IS AN APPLICATION OF FOURTH-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
C...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
C...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
C...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
C...  GAZETTE, PP. 19-27, 1941, N = 4, M = 1, P = 0, 1, 2, 3, 4.  THE
C...  IMPLEMENTATION IS THE **FIVE-POINT BIASED UPWIND FORMULA** OF
C...  M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE
C...  ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69,
C...  AUGUST, 1978
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION      DX,  R4FDX,      U,     UX,      V,     XL,
     1                      XU, DFLOAT
C...
      DIMENSION U(N),UX(N)
C...
C...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
C...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
C...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      DX=(XU-XL)/DFLOAT(N-1)
      R4FDX=1.D+00/(12.D+00*DX)
      IF(V.LT.0.D+00)GO TO 10
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
      UX(  1)=R4FDX*
     1( -25.D+00     *U(  1)
     2  +48.D+00     *U(  2)
     3  -36.D+00     *U(  3)
     4  +16.D+00     *U(  4)
     5   -3.D+00     *U(  5))
      UX(  2)=R4FDX*
     1(  -3.D+00     *U(  1)
     2  -10.D+00     *U(  2)
     3  +18.D+00     *U(  3)
     4   -6.D+00     *U(  4)
     5   +1.D+00     *U(  5))
      UX(  3)=R4FDX*
     1(  +1.D+00     *U(  1)
     2   -8.D+00     *U(  2)
     3   +0.D+00     *U(  3)
     4   +8.D+00     *U(  4)
     5   -1.D+00     *U(  5))
      NM1=N-1
      DO 1 I=4,NM1
      UX(  I)=R4FDX*
     1(  -1.D+00     *U(I-3)
     2   +6.D+00     *U(I-2)
     3  -18.D+00     *U(I-1)
     4  +10.D+00     *U(I  )
     5   +3.D+00     *U(I+1))
1     CONTINUE
      UX(  N)=R4FDX*
     1(  +3.D+00     *U(N-4)
     2  -16.D+00     *U(N-3)
     3  +36.D+00     *U(N-2)
     4  -48.D+00     *U(N-1)
     5  +25.D+00     *U(N  ))
      RETURN
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
10    UX(  1)=R4FDX*
     1( -25.D+00     *U(  1)
     2  +48.D+00     *U(  2)
     3  -36.D+00     *U(  3)
     4  +16.D+00     *U(  4)
     5   -3.D+00     *U(  5))
      NM3=N-3
      DO 2 I=2,NM3
      UX(  I)=R4FDX*
     1(  -3.D+00     *U(I-1)
     2  -10.D+00     *U(I  )
     3  +18.D+00     *U(I+1)
     4   -6.D+00     *U(I+2)
     5   +1.D+00     *U(I+3))
2     CONTINUE
      UX(N-2)=R4FDX*
     1(  +1.D+00     *U(N-4)
     2   -8.D+00     *U(N-3)
     3   +0.D+00     *U(N-2)
     4   +8.D+00     *U(N-1)
     5   -1.D+00     *U(N  ))
      UX(N-1)=R4FDX*
     1(  -1.D+00     *U(N-4)
     2   +6.D+00     *U(N-3)
     3  -18.D+00     *U(N-2)
     4  +10.D+00     *U(N-1)
     5   +3.D+00     *U(N  ))
      UX(  N)=R4FDX*
     1(  +3.D+00     *U(N-4)
     2  -16.D+00     *U(N-3)
     3  +36.D+00     *U(N-2)
     4  -48.D+00     *U(N-1)
     5  +25.D+00     *U(N  ))
      RETURN
      END
*DECK DSS024
      SUBROUTINE DSS024(XL,XU,N,U,UX,V)
C...
C...  SUBROUTINE DSS024 IS AN APPLICATION OF FOURTH-ORDER DIRECTIONAL
C...  DIFFERENCING IN THE NUMERICAL METHOD OF LINES.  IT IS INTENDED
C...  SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY
C...  FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS-
C...  CUSSED IN SUBROUTINE DSS012.  THE COEFFICIENTS OF THE FINITE
C...  DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W.
C...  G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL
C...  GAZETTE, PP. 19-27, 1941, N = 4, M = 1, P = 0, 1, 2, 3, 4.  THE
C...  IMPLEMENTATION IS THE **FIVE-POINT BIASED UPWIND FORMULA** OF
C...  M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE
C...  ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69,
C...  AUGUST, 1978.
C...
C...  DSS024 IS AN EXTENSION OF DSS020.  THESE ROUTINES HAVE THE 
C...  FOLLOWING APPROXIMATIONS:
C...
C...     DSS020 - FIVE POINT APPROXIMATIONS OVER THE ENTIRE SPATIAL
C...              GRID.  AT THE INTERIOR POINTS, THE FIVE POINT, BIASED
C...              UPWIND APPROXMTAION IS USED.  
C...
C...     DSS024 - A COMBINATION OF TWO POINT, UPWIND AND FIVE POINT, 
C...              BIASED UPWIND APPROXIMATIONS IS USED.  THE FIVE POINT
C...              APPROXIMATION IS USED AT THE INTERIOR POINTS TO
C...              PROVIDE RESOLUTION OF MOVING FRONTS (I.E., TO MINIMIZE
C...              NUMERICAL DIFFUSION).  THE TWO POINT APPROXIMATION 
C...              IS USED NEAR THE BOUNDARIES TO REDUCE NUMERICAL 
C...              OSCILLATIONS FROM THE FIVE POINT APPROXIMATION.
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(N),UX(N)
C...
C...  NUMBER OF GRID POINTS AT THE LEFT AND RIGHT ENDS WITH TWO POINT
C...  UPWIND APPROXIMATIONS
      NL=5
      NU=5
C...
C...  COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION
C...  CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER-
C...  ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT).
      DX=(XU-XL)/DFLOAT(N-1)
      R4FDX=1.D+00/(12.D+00*DX)
      IF(V.GE.0.D+00)THEN
C...
C...     (1)  FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V
C...
C...  TWO POINT UPWIND APPROXIMATIONS AT LEFT END (FOR NL GRID POINTS)
      UX(  1)=(1.0D0/DX)*
     1(  -1.D+00     *U(1)
     2   +1.D+00     *U(2))
      DO 2 I=2,NL
      UX(  I)=(1.0D0/DX)*
     1(  -1.D+00     *U(I-1)
     2   +1.D+00     *U(I  ))
2     CONTINUE
C...
C...  FIVE POINT BIASED UPWIND APPROXIMATIONS FOR INTERIOR GRID POINTS 
C...  NL+1 TO N-NU-1
      DO 1 I=NL+1,N-NU-1
      UX(  I)=R4FDX*
     1(  -1.D+00     *U(I-3)
     2   +6.D+00     *U(I-2)
     3  -18.D+00     *U(I-1)
     4  +10.D+00     *U(I  )
     5   +3.D+00     *U(I+1))
1     CONTINUE
C...
C...  TWO POINT UPWIND APPROXIMATIONS AT RIGHT END (FOR NU GRID POINTS)
      DO 3 I=N-NU,N
      UX(I  )=(1.0D0/DX)*
     1(  -1.D+00     *U(I-1)
     2   +1.D+00     *U(I  ))
3     CONTINUE
C...
C...     (2)  FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V
C...
      ELSE IF(V.LT.0.D+00)THEN
C...
C...  TWO POINT UPWIND APPROXIMATIONS AT LEFT END (FOR NL GRID POINTS)
      DO 12 I=1,NL
      UX(  I)=(1.0D0/DX)*
     1(  -1.D+00     *U(I  )
     2   +1.D+00     *U(I+1))
12    CONTINUE
C...
C...  FIVE POINT BIASED UPWIND APPROXIMATIONS FOR INTERIOR GRID POINTS 
C...  NL+1 TO N-NU-1
      DO 11 I=NL+1,N-NU-1
      UX(  I)=R4FDX*
     1(  -3.D+00     *U(I-1)
     2  -10.D+00     *U(I  )
     3  +18.D+00     *U(I+1)
     4   -6.D+00     *U(I+2)
     5   +1.D+00     *U(I+3))
11    CONTINUE
C...
C...  TWO POINT UPWIND APPROXIMATIONS AT RIGHT END (FOR NU GRID POINTS)
      DO 13 I=N-NU,N-1
      UX(  I)=(1.0D0/DX)*
     1(  -1.D+00     *U(I  )
     2   +1.D+00     *U(I+1))
13    CONTINUE
      UX(  N)=(1.0D0/DX)*
     1(  -1.D+00     *U(N-1)
     2   +1.D+00     *U(N  ))
      END IF
      RETURN
      END
*DECK DSS026
      SUBROUTINE DSS026(XL,XU,N1,N2,ND,U2D,UX2D,V)
C...
C...  SUBROUTINE DSS026 COMPUTES A PARTIAL DERIVATIVE OVER A TWO-
C...  DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE-
C...  POINT BIASED UPWIND APPROXIMATIONS.  IT IS INTENDED PRIMARILY
C...  FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION
C...  OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN TWO DIMENSIONS.  IT
C...  DIFFERS FROM DSS034 ONLY IN THAT DSS024 IS CALLED IN PLACE OF
C...  DSS020.  IT THEREFORE SUPERSEDES DSS034 (IN THE SENSE THAT DSS024
C...  SUPERSEDES DSS020).
C...
C...  ARGUMENT LIST
C...
C...     XL        LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     XU        UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     N1        NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     N2        NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     ND        NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE
C...               PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     U2D       TWO-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT VARI-
C...               ABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (INPUT)
C...
C...     UX2D      TWO-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI-
C...               VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (OUTPUT)
C...
C...     V         VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED
C...               OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE
C...               PARTIAL DERIVATIVE.  V EQ 0 CALLS THE FIVE-POINT
C...               CENTERED APPROXIMATION.  V NE 0 CALLS THE FIVE-POINT
C...               BIASED UPWIND APPROXIMATION (INPUT)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION   UX1D,   UX2D,    U1D,    U2D,      V,     XL,
     1                     XU
C...
C...  THE FOLLOWING TWO-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U2D) AND ITS PARTIAL DERIVATIVE (UX2D)
      DIMENSION   U2D(N1,N2), UX2D(N1,N2)
C...
C...  THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D).  IN EACH
C...  CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE
C...  OTHER INDEPENDENT VARIABLE VARIES OVER ITS TOTAL INTERVAL.
C...  THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE
C...  ONE-DIMENSIONAL ROUTINES DSS004 AND DSS024.
C...
C...  NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE-
C...  FORE HAVE TO BE INCREASED IN SIZE.  HOWEVER, WITH A SIZE
C...  OF 101, THE TWO-DIMENSIONAL PROBLEM COULD HAVE A GRID OF
C...  101 X 101 POINTS, THEREBY GENERATING AN APPROXIMATING ODE
C...  SYSTEM WITH A MULTIPLE OF 101 X 101 EQUATIONS, DEPENDING ON
C...  THE NUMBER OF SIMULTANEOUS PDES.  THIS IS A VERY LARGE ODE
C...  PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING
C...  IS CONSIDERED ADEQUATE FOR MOST PROBLEMS.
      PARAMETER (NSIZE=101)
      DIMENSION U1D(NSIZE), UX1D(NSIZE)
C...
C...  GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED
C...  WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE
      IF(ND.EQ.2)GO TO 2
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N1 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 10, 11 AND 12
      DO 10 J=1,N2
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS024 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 11 I=1,N1
      U1D(I)=U2D(I,J)
11    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS024(XL,XU,N1,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UX2D
      DO 12 I=1,N1
      UX2D(I,J)=UX1D(I)
12    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE SECOND INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE SECOND INDEPENDENT VARIABLE
10    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D
      RETURN
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N2 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 20 AND 21
2     DO 20 I=1,N1
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS024 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 21 J=1,N2
      U1D(J)=U2D(I,J)
21    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS024(XL,XU,N2,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UX2D
      DO 22 J=1,N2
      UX2D(I,J)=UX1D(J)
22    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE FIRST INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE FIRST INDEPENDENT VARIABLE
20    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D
      RETURN
      END
*DECK DSS028
      SUBROUTINE DSS028(XL,XU,N1,N2,ND,U2D,UX2D,UXX2D,NL,NU)
C...
C...  SUBROUTINE DSS028 COMPUTES A SECOND ORDER PARTIAL DERIVATIVE
C...  OVER A TWO-DIMENSIONAL DOMAIN USING FIVE-POINT CENTERED
C...  APPROXIMATIONS.  IT IS INTENDED PRIMARILY FOR THE NUMERICAL
C...  INTEGRATION OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN TWO
C...  DIMENSIONS.
C...
C...  WRITTEN BY JOHN SALMENTO, DEC, 1992, WITH MINOR EDITING BY 
C...  W. E. SCHIESSER.
C...
C...  DSS028 CALLS DSS044 TO COMPUTE THE SECOND DERIVATIVE (WITHIN
C...  A SERIES OF LOOPS TO COVER A TWO-DIIMENSIONAL DOMAIN, RATHER
C...  THAN THE ONE-DIMENSIONAL DOMAIN COVERED BY DSS044).
C...
C...  ARGUMENT LIST
C...
C...     XL        LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     XU        UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     N1        NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     N2        NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     ND        NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE
C...               PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     U2D       TWO-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT VARI-
C...               ABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (INPUT)
C...
C...     UX2D      TWO-DIMENSIONAL ARRAY CONTAINING THE FIRST ORDER
C...               PARTIAL DERIVATIVE OF THE DEPENDENT VARIABLE WITH 
C...               RESPECT TO INDEPENDENT VARIABLE ND (INPUT)
C...
C...     UXX2D     TWO-DIMENSIONAL ARRAY CONTAINING THE SECOND ORDER
C...               PARTIAL DERIVATIVE OF THE DEPENDENT VARIABLE WITH 
C...               RESPECT TO INDEPENDENT VARIABLE ND (OUTPUT)
C...
C...     NL        INTEGER ARRAY INDEX FOR THE TYPE OF BOUNDARY CONDITION
C...               AT X = XL (INPUT)  THE ALLOWABLE VALUES ARE
C...
C...                  1 - DIRICHLET BOUNDARY CONDITION AT X = XL
C...
C...                      UX2D(*,1) IS NOT USED IF ND = 1
C...
C...                      UX2D(1,*) IS NOT USED IF ND = 2
C...
C...                  2 - NEUMANN BOUNDARY CONDITION AT X = XL
C...
C...                      UX2D(*,1) IS USED IF ND = 1
C...
C...                      UX2D(1,*) IS USED IF ND = 2
C...
C...     NU        INTEGER ARRAY INDEX FOR THE TYPE OF BOUNDARY CONDITION
C...               AT X = XU (INPUT)  THE ALLOWABLE VALUES ARE
C...
C...                  1 - DIRICHLET BOUNDARY CONDITION AT X = XU
C...
C...                      UX2D(*,N2) IS NOT USED IF ND = 1
C...
C...                      UX2D(N1,*) IS NOT USED IF ND = 2
C...
C...                  2 - NEUMANN BOUNDARY CONDITION AT X = XU
C...
C...                      UX2D(*,N2) IS USED IF ND = 1
C...
C...                      UX2D(N1,*) IS USED IF ND = 2
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  ACCESS THE INPUT/OUTPUT UNIT NUMBERS FOR PRINTING ERROR
C...  MESSAGES
      COMMON/IO/NI, NO
C...
C...  THE FOLLOWING TWO-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U2D) AND ITS PARTIAL DERIVATIVES (UX2D, UXX2D)
      DIMENSION U2D(N1,N2), UX2D(N1,N2), UXX2D(N1,N2)
C...
C...  DIMENSION THE INTEGER ARRAYS CONTAINING THE BOUNDARY CON-
C...  DITION INDICES
      DIMENSION NL(*), NU(*)
C...
C...  THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U1D) AND ITS PARTIAL DERIVATIVES (UX1D, UXX1D)).
C...  IN EACH CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT 
C...  AND THE OTHER INDEPENDENT VARIABLE VARIES OVER ITS TOTAL 
C...  INTERVAL.  THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN 
C...  CALLING THE ONE-DIMENSIONAL ROUTINE DSS044.  NOTE THAT THE 
C...  ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THEREFORE HAVE TO 
C...  BE INCREASED IN SIZE.
      PARAMETER (NSIZE=101)
      DIMENSION U1D(NSIZE), UX1D(NSIZE), UXX1D(NSIZE)
C...
C...  GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED
C...  WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE
      IF(ND.EQ.2)GO TO 2
C...
C...  ******************************************************************
C...
C...  TEST IF N1 IS GREATER THAN THE DIMENSION OF THE TEMPORARY ARRAYS
C...  U1D, UX1D AND UXX1D.  IF IT IS, WRITE A WARNING MESSAGE TO UNIT
C...  NO AND STOP EXECUTION SINCE SUBROUTINE DSS028 WILL HAVE AN ARRAY
C...  OVERFLOW
      IF(N1.GT.NSIZE)THEN
         WRITE(*,31)NSIZE
31       FORMAT(' Size of first dimension (N1) in arrays passed to ',/,
     +   'DSS028 is larger than the maximum size ',I3,'.  Increase ',/,
     +   'NSIZE in DSS028 and recompile this subroutine.')
      STOP
      END IF
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N1 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 10, 11 AND 12
      DO 10 J=1,N2
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAYS 
C...  U2D AND UX2D TO THE ONE-DIMENSIONAL ARRAYS U1D AND UX1D SO 
C...  THAT SUBROUTINE DSS044 CAN BE USED TO CALCULATE THE SECOND
C...  ORDER PARTIAL DERIVATIVE
      DO 11 I=1,N1
          U1D(I)= U2D(I,J)
         UX1D(I)=UX2D(I,J)
11    CONTINUE
C...
C...  COMPUTE THE SECOND DERIVATIVE UXX1D
      CALL DSS044(XL,XU,N1,U1D,UX1D,UXX1D,NL(J),NU(J))
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UXX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UXX2D
      DO 12 I=1,N1
         UXX2D(I,J)=UXX1D(I)
12    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE SECOND INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE SECOND INDEPENDENT VARIABLE
10    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UXX2D
      RETURN
C...
C...  ******************************************************************
C...
C...  TEST IF N2 IS GREATER THAN THE DIMENSION OF THE TEMPORARY ARRAYS
C...  U1D, UX1D AND UXX1D.  IF IT IS, WRITE A WARNING MESSAGE TO UNIT
C...  NO AND STOP EXECUTION SINCE SUBROUTINE DSS028 WILL HAVE AN ARRAY
C...  OVERFLOW
      IF(N2.GT.NSIZE)THEN
         WRITE(*,32)NSIZE
32       FORMAT(' Size of second dimension (N2) in arrays passed to ',/,
     +   'DSS028 is larger than the maximum size ',I3,'.  Increase ' ,/,
     +   'NSIZE in DSS028 and recompile this subroutine.')
      STOP
      END IF
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N2 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 20 AND 21
2     DO 20 I=1,N1
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAYS 
C...  U2D AND UX2D TO THE ONE-DIMENSIONAL ARRAYS U1D AND UX1D SO 
C...  THAT SUBROUTINE DSS044 CAN BE USED TO CALCULATE THE SECOND
C...  ORDER PARTIAL DERIVATIVE
      DO 21 J=1,N2
          U1D(J)= U2D(I,J)
         UX1D(J)=UX2D(I,J)
21    CONTINUE
C...
C...  COMPUTE THE SECOND DERIVATIVE UXX1D
      CALL DSS044(XL,XU,N2,U1D,UX1D,UXX1D,NL(I),NU(I))
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UXX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UXX2D
      DO 22 J=1,N2
         UXX2D(I,J)=UXX1D(J)
22    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE FIRST INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE FIRST INDEPENDENT VARIABLE
20    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UXX2D
      RETURN
      END
*DECK DSS032
      SUBROUTINE DSS032(NCALL,NPTS,IGRID,N,X,U,UX,C)
C...
C...  SUBROUTINE DSS032 COMPUTES A NUMERICAL DERIVATIVE BY FIVE-POINT
C...  CENTERED AND NONCENTERED APPROXIMATIONS BASED ON THE FOURTH-ORDER
C...  LAGRANGE INTERPOLATON POLYNOMIAL.  THE MATHEMATICAL DETAILS OF
C...  THE APPROXIMATIONS ARE GIVEN IN THE SUBSEQUENT COMMENTS, AND IN
C...  THE COMMENTS IN SUBROUTINE DSS32A, THE CORE DIFFERENTIATOR CALLED
C...  BY DSSO32.
C...
C...  ARGUMENT LIST
C...
C...     NCALL   INTEGER INDEX FOR SPECIFYING THE CALCULATION OF THE
C...             WEIGHTING COEFFICIENTS IN THE DIFFERENTIATION FORMULAS
C...             DURING THE FIRST CALL TO SUBROUTINE DSS032.  SPECIFI-
C...             CALLY,
C...
C...                   NCALL = 1 -  THE FIRST CALL TO SUBROUTINE DSS032
C...                                FOR WHICH THE WEIGHTING COEFFICIENTS
C...                                ARE CALCULATED
C...
C...                   NCALL NE 1 - A CALL TO SUBROUTINE DSS032 AFTER
C...                                THE FIRST CALL SO THAT THE WEIGHTING
C...                                COEFFICIENTS ARE NOT RECALCULATED
C...
C...             (INPUT)
C...
C...     NPTS    NUMBER OF POINTS USED IN THE NUMERICAL DIFFERENTIATION
C...             FORMULAS.  PRESENTLY, ONLY FIVE-POINT FORMULAS ARE USED
C...             IN SUBROUTINE DSS032 SO NPTS = 5 (INPUT)
C...
C...     IGRID   ONE-DIMENSIONAL ARRAY CONTAINING THE IDENTIFIERS FOR
C...             THE FIVE BASIC GRIDS DESCRIBED ABOVE (INPUT)
C...
C...     N       NUMBER OF GRID POINTS OVER WHICH THE INDEPENDENT
C...             VARIABLE, X, AND DEPENDENT VARIABLE, U, ARE DEFINED,
C...             INCLUDING THE END POINTS (INPUT)
C...
C...     X       ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE
C...             INDEPENDENT VARIABLE (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE
C...             DEPENDENT VARIABLE (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE
C...             NUMERICAL DERIVATIVE OF U WITH RESPECT TO X (OUTPUT)
C...
C...     C       TWO-DIMENSIONAL ARRAY CONTAINING THE WEIGHTING COEFFI-
C...             CIENTS IN THE DIFFERENTIATION FORMULAS (OUTPUT)
C...
C...  SUBROUTINE DSS032 IS A NUMERICAL DIFFERENTIATOR BASED ON THE FIVE-
C...  POINT LAGRANGIAN INTERPOLATION POLYNOMIAL.  A NUMERICAL DERIVATIVE
C...  UX IS COMPUTED BASED ON ONE OF THE FOLLOWING FIVE-POINT GRIDS
C...
C...     (1)      UX
C...               I    I+1   I+2   I+3   I+4
C...
C...     (2)            UX
C...              I-1    I    I+1   I+2   I+3
C...
C...     (3)                  UX
C...              I-2   I-1    I    I+1   I+2
C...
C...     (4)                        UX
C...              I-3   I-2   I-1    I    I+1
C...
C...     (5)                              UX
C...              I-4   I-3   I-2   I-1    I
C...
C...  FOR GRID (1), FOUR POINTS ARE USED TO THE RIGHT OF THE POINT WHERE
C...  THE NUMERICAL DERIVATIVE, UX, IS COMPUTED.  FOR GRID (2), ONE
C...  POINT TO THE LEFT AND THREE POINTS TO THE RIGHT OF THE POINT WHERE
C...  THE NUMERICAL DERIVATIVE IS COMPUTED ARE USED, ETC.  THUS, GRID
C...  (3) IS A FIVE-POINT CENTERED APPROXIMATION, GRIDS (2) AND (4) ARE
C...  FIVE-POINT BIASED UPWIND APPROXIMATIONS, AND GRIDS (1) AND (5) ARE
C...  FIVE-POINT UPWIND APPROXIMATIONS.
C...
C...  THE USER MAY SELECT ANY COMBINATION OF THESE GRIDS OVER A DOMAIN
C...  XL LE X LE XU CONSISTING OF A TOTAL OF N GRID POINTS.  NORMALLY,
C...  GRIDS (1) AND (5) ARE USED AT THE LEFT AND RIGHT BOUNDARIES, X =
C...  XL AND X = XU, RESPECTIVELY.  GRIDS (2) AND (4) ARE USED AT ONE
C...  GRID POINT TO THE RIGHT AND LEFT OF THE LEFT AND RIGHT BOUNDARIES,
C...  RESPECTIVELY.
C...
C...  GRID (3) CAN BE USED IN THE INTERIOR OF THE X DOMAIN IN PROBLEMS
C...  FOR WHICH A CENTERED APPROXIMATION IS APPROPRIATE, E.G., A PARA-
C...  BOLIC PROBLEM IN PARTIAL DIFFERENTIAL EQUATIONS (PDES).
C...
C...  GRID (4) CAN BE USED THROUGHOUT THE INTERIOR OF THE X DOMAIN FOR
C...  FIRST-ORDER HYPERBOLIC PDES WITH A POSITIVE VELOCITY, V, IN THE
C...  DIFFERENTIAL GROUP UT + V*UX, I.E., BIASED UPWIND APPROXIMATION OF
C...  THE SPATIAL DERIVATIVE, UX, WITH GREATER WEIGHTING TO THE LEFT
C...  (UPWIND) OF THE POINT WHERE UX IS CALCULATED.  CONVERSELY, GRID
C...  (2) CAN BE USED THROUGHOUT THE INTERIOR OF THE X DOMAIN FOR FIRST-
C...  ORDER HYPERBOLIC PDES WITH A NEGATIVE VELOCITY, I.E., BIASED UP-
C...  WIND APPROXIMATION WITH GREATER WEIGHTING TO THE RIGHT (UPWIND) OF
C...  THE POINT WHERE UX IS CALCULATED.
C...
C...  THE SELECTION OF THE FIVE BASIC GRIDS IS MADE BY THE USER THROUGH
C...  GRID IDENTIFIER ARRAY IGRID, THE THIRD ARGUMENT OF SUBROUTINE
C...  DSS032.  FOR EXAMPLE, IF IGRID(1) = 1, IGRID(2) = 2, IGRID(I) = 3,
C...  IGRID(N-1) = 4, IGRID(N) = 5, THE FIRST AND SECOND GRID POINTS
C...  WILL RETURN THE DERIVATIVE UX ACCORDING TO GRIDS (1) AND (2)
C...  ABOVE.  THE (N-1)ST AND NTH GRID POINTS WILL RETURN UX ACCORDING
C...  TO GRIDS (4) AND (5) ABOVE.  THE INTERIOR GRID POINTS, WITH SUB-
C...  SCRIPT I, 3 LE I LE N-2, WILL RETURN UX ACCORDING TO GRID (3)
C...  ABOVE.
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION      C,     C5,      U,     UX,      X,     X5
C...
C...  VARIABLE DIMENSION THE ARRAYS FOR THE ARGUMENTS LISTED ABOVE
      DIMENSION IGRID(N), X(N), U(N), UX(N), C(NPTS,N)
C...
C...  ABSOLUTE DIMENSION THE ARRAYS USED BY SUBROUTINE DSS32A WHICH IS
C...  CALLED BY DSS032
      DIMENSION X5(5), C5(5)
C...
C...  IF THIS IS THE FIRST CALL TO SUBROUTINE DSS032, GO TO STATEMENT
C...  17 TO CALCULATE THE COEFFICIENTS IN THE DIFFERENTIATION FORMULAS
C...  VIA SUBROUTINE DSS32A
      IF(NCALL.EQ.1)GO TO 17
C...
C...  THE NUMERICAL DERIVATIVES ARE CALCULATED ON THE N GRID POINTS
C...  USING THE WEIGHTING COEFFICIENTS CALCULATED PREVIOUSLY BY SUB-
C...  ROUTINE DSS32A
7     DO 6 I=1,N
      ITYPE=IGRID(I)
      GO TO(1,2,3,4,5),ITYPE
C...
C...  GRID TYPE 1
1     UX(I)=C(1,I)*U(I)
     1     +C(2,I)*U(I+1)
     2     +C(3,I)*U(I+2)
     3     +C(4,I)*U(I+3)
     4     +C(5,I)*U(I+4)
      GO TO 6
C...
C...  GRID TYPE 2
2     UX(I)=C(1,I)*U(I-1)
     1     +C(2,I)*U(I)
     2     +C(3,I)*U(I+1)
     3     +C(4,I)*U(I+2)
     4     +C(5,I)*U(I+3)
      GO TO 6
C...
C...  GRID TYPE 3
3     UX(I)=C(1,I)*U(I-2)
     1     +C(2,I)*U(I-1)
     2     +C(3,I)*U(I)
     3     +C(4,I)*U(I+1)
     4     +C(5,I)*U(I+2)
      GO TO 6
C...
C...  GRID TYPE 4
4     UX(I)=C(1,I)*U(I-3)
     1     +C(2,I)*U(I-2)
     2     +C(3,I)*U(I-1)
     3     +C(4,I)*U(I)
     4     +C(5,I)*U(I+1)
      GO TO 6
C...
C...  GRID TYPE 5
5     UX(I)=C(1,I)*U(I-4)
     1     +C(2,I)*U(I-3)
     2     +C(3,I)*U(I-2)
     3     +C(4,I)*U(I-1)
     4     +C(5,I)*U(I)
6     CONTINUE
      RETURN
C...
C...  DURING THE FIRST CALL TO SUBROUTINE DSS032, THE COEFFICIENTS IN
C...  THE DIFFERENTIATION FORMULAS, STORED IN ARRAY C, ARE CALCULATED
C...  BY A CALL TO SUBROUTINE DSS32A
17    DO 16 I=1,N
C...
C...  DEFINE THE VALUES OF X USED BY SUBROUTINE DSS32A TO CALCULATE
C...  THE WEIGHTING COEFFICIENTS IN THE DIFFERENTIATION FORMULAS
      ITYPE=IGRID(I)
      GO TO(11,12,13,14,15),ITYPE
C...
C...  GRID TYPE 1
11    X5(1)=X(I)
      X5(2)=X(I+1)
      X5(3)=X(I+2)
      X5(4)=X(I+3)
      X5(5)=X(I+4)
      GO TO 19
C...
C...  GRID TYPE 2
12    X5(1)=X(I-1)
      X5(2)=X(I)
      X5(3)=X(I+1)
      X5(4)=X(I+2)
      X5(5)=X(I+3)
      GO TO 19
C...
C...  GRID TYPE 3
13    X5(1)=X(I-2)
      X5(2)=X(I-1)
      X5(3)=X(I)
      X5(4)=X(I+1)
      X5(5)=X(I+2)
      GO TO 19
C...
C...  GRID TYPE 4
14    X5(1)=X(I-3)
      X5(2)=X(I-2)
      X5(3)=X(I-1)
      X5(4)=X(I)
      X5(5)=X(I+1)
      GO TO 19
C...
C...  GRID TYPE 5
15    X5(1)=X(I-4)
      X5(2)=X(I-3)
      X5(3)=X(I-2)
      X5(4)=X(I-1)
      X5(5)=X(I)
C...
C...  CALL SUBROUTINE DSS32A TO CALCULATE THE WEIGHTING COEFFICIENTS IN
C...  DIFFERENTIATION FORMULA FOR GRID POINT I
19    CALL DSS32A(ITYPE,X5,C5)
C...
C...  STORE THE FIVE WEIGHTING COEFFICIENTS IN ARRAY C5 IN THE TWO-
C...  DIMENSIONAL ARRAY C
      DO 18 IC=1,5
      C(IC,I)=C5(IC)
18    CONTINUE
16    CONTINUE
C...
C...  ALL OF THE WEIGHTING COEFFICIENTS FOR THE N GRID POINTS HAVE BEEN
C...  CALCULATED AND STORED IN ARRAY C SO RETURN TO DO LOOP 6 TO CALCU-
C...  LATE THE NUMERICAL DERIVATIVES
      GO TO 7
      END
      SUBROUTINE DSS32A(I,X,C)
C...
C...  SUBROUTINE DSS32A IS THE CORE DIFFERENTIATOR CALLED BY SUBROUTINE
C...  DSS032.  IT BASICALLY COMPUTES THE WEIGHTING COEFFCIENTS IN THE
C...  FIVE-POINT CENTERED AND NONCENTERED DIFFERENTIATION FORMULAS FOR
C...  USE IN SUBROUTINE DSS032.
C...
C...  ARGUMENT LIST
C...
C...     I       TYPE OF GRID AS DISCUSSED AT THE BEGINNING OF SUBROU-
C...             TINE DSS032 WITH POSSIBLE VALUES 1, 2, 3, 4, 5 (INPUT)
C...
C...     X       ONE-DIMENSIONAL ARRAY CONTAINING THE FIVE VALUES OF THE
C...             INDEPENDENT VARIABLE, X, USED IN THE CALCULATION OF THE
C...             NUMERICAL DERIVATIVE (INPUT)
C...
C...     C       ONE-DIMENSIONAL ARRAY CONTAINING THE FIVE COEFFICIENTS
C...             OF THE DIFFERENTIATION FORMULA FOR GRID TYPE I (OUTPUT)
C...
C...  SUBROUTINE DSS32A COMPUTES THE WEIGHTING COEFFICIENTS IN THE
C...  DIFFERENTIATION FORMULAS USED IN SUBROUTINE DSS032.  DSS032
C...  CALLS DSS32A JUST ONCE TO CALCULATE THE WEIGHTING COEFFICIENTS
C...  AFTER THE VARIABLE GRID IS DEFINED SINCE THE WEIGHTING CO-
C...  EFFICIENTS ARE CONSTANTS FOR A GIVEN VARIABLE GRID.
C...
C...  THE GENERAL FIVE-POINT LAGRANGE INTERPOLATION POLYNOMIAL SERVES
C...  AS THE MATHEMATICAL BASIS FOR THIS SUBROUTINE
C...
C...
C...            (X  - X2)(X  - X3)(X  - X4)(X  - X5)
C...     U(X) = .................................... U(X1)
C...            (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...            (X  - X1)(X  - X3)(X  - X4)(X  - X5)
C...          + .................................... U(X2)
C...            (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...            (X  - X1)(X  - X2)(X  - X4)(X  - X5)
C...          + .................................... U(X3)           (1)
C...            (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...            (X  - X1)(X  - X2)(X  - X3)(X  - X5)
C...          + .................................... U(X4)
C...            (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...            (X  - X1)(X  - X2)(X  - X3)(X  - X4)
C...          + .................................... U(X5)
C...            (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...  DIFFERENTIATION OF EQUATION (1) GIVES THE FIVE-POINT DIFFERENTIA-
C...  TION FORMULAS USED IN SUBROUTINE DSS032 FOR THE CALCULATION OF
C...  NUMERICAL DERIVATIVES OVER A VARIABLE GRID
C...
C...
C...     UX(X) = ((X - X2)(X - X3)(X - X4) +
C...
C...              (X - X2)(X - X3)(X - X5) +
C...
C...              (X - X2)(X - X4)(X - X5) +
C...
C...              (X - X3)(X - X4)(X - X5))
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...           + ((X - X1)(X - X3)(X - X4) +
C...
C...              (X - X1)(X - X3)(X - X5) +
C...
C...              (X - X1)(X - X4)(X - X5) +
C...
C...              (X - X3)(X - X4)(X - X5))
C...         .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...           + ((X - X1)(X - X2)(X - X4) +
C...
C...              (X - X1)(X - X2)(X - X5) +
C...
C...              (X - X1)(X - X4)(X - X5) +
C...
C...              (X - X2)(X - X4)(X - X5))
C...         .................................... U(X3)              (2)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...           + ((X - X1)(X - X2)(X - X3) +
C...
C...              (X - X1)(X - X2)(X - X5) +
C...
C...              (X - X1)(X - X3)(X - X5) +
C...
C...              (X - X2)(X - X3)(X - X5))
C...         .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...           + ((X - X1)(X - X2)(X - X3) +
C...
C...              (X - X1)(X - X2)(X - X4) +
C...
C...              (X - X1)(X - X3)(X - X4) +
C...
C...              (X - X2)(X - X3)(X - X4))
C...         .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...  IF X = X1 IN EQUATION (2), THE FORMULA FOR THE FIRST DERIVATIVE
C...  AT X = X1, DENOTED AS UX(X1), RESULTS
C...
C...
C...    UX(X1) = ((X1 - X2)(X1 - X3)(X1 - X4) +
C...
C...              (X1 - X2)(X1 - X3)(X1 - X5) +
C...
C...              (X1 - X2)(X1 - X4)(X1 - X5) +
C...
C...              (X1 - X3)(X1 - X4)(X1 - X5))
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...             ((X1 - X3)(X1 - X4)(X1 - X5))
C...       + .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...             ((X1 - X2)(X1 - X4)(X1 - X5))
C...       + .................................... U(X3)              (3)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...             ((X1 - X2)(X1 - X3)(X1 - X5))
C...       + .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...             ((X1 - X2)(X1 - X3)(X1 - X4))
C...       + .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...  SIMILARLY, THE SUBSTITUTIONS X = X2, X = X3, X = X4, X = X5 IN
C...  EQUATION (2) GIVES RESPECTIVELY THE DIFFERENTIATION FORMULAS FOR
C...  UX(X2), UX(X3), UX(X4), UX(X5)
C...
C...
C...    UX(X2) = ((X2 - X3)(X2 - X4)(X2 - X5))
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...           + ((X2 - X1)(X2 - X3)(X2 - X4) +
C...
C...              (X2 - X1)(X2 - X3)(X2 - X5) +
C...
C...              (X2 - X1)(X2 - X4)(X2 - X5) +
C...
C...              (X2 - X3)(X2 - X4)(X2 - X5))
C...         .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...             ((X2 - X1)(X2 - X4)(X2 - X5))
C...       + .................................... U(X3)              (4)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...             ((X2 - X1)(X2 - X3)(X2 - X5))
C...       + .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...             ((X2 - X1)(X2 - X3)(X2 - X4))
C...       + .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...
C...    UX(X3) = ((X3 - X2)(X3 - X4)(X3 - X5))
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...             ((X3 - X1)(X3 - X4)(X3 - X5))
C...       + .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...           + ((X3 - X1)(X3 - X2)(X3 - X4) +
C...
C...              (X3 - X1)(X3 - X2)(X3 - X5) +
C...
C...              (X3 - X1)(X3 - X4)(X3 - X5) +
C...
C...              (X3 - X2)(X3 - X4)(X3 - X5))
C...         .................................... U(X3)              (5)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...             ((X3 - X1)(X3 - X2)(X3 - X5))
C...       + .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...             ((X3 - X1)(X3 - X2)(X3 - X4))
C...       + .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...
C...    UX(X4) = ((X4 - X2)(X4 - X3)(X4 - X5))
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...             ((X4 - X1)(X4 - X3)(X4 - X5))
C...       + .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...             ((X4 - X1)(X4 - X2)(X4 - X5))
C...       + .................................... U(X3)              (6)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...           + ((X4 - X1)(X4 - X2)(X4 - X3) +
C...
C...              (X4 - X1)(X4 - X2)(X4 - X5) +
C...
C...              (X4 - X1)(X4 - X3)(X4 - X5) +
C...
C...              (X4 - X2)(X4 - X3)(X4 - X5))
C...         .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...             ((X4 - X1)(X4 - X2)(X4 - X3))
C...       + .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...
C...    UX(X5) = ((X5 - X2)(X5 - X3)(X5 - X4) +
C...         .................................... U(X1)
C...         (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5)
C...
C...
C...             ((X5 - X1)(X5 - X3)(X5 - X4))
C...       + .................................... U(X2)
C...         (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5)
C...
C...
C...             ((X5 - X1)(X5 - X2)(X5 - X4))
C...       + .................................... U(X3)              (7)
C...         (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5)
C...
C...
C...             ((X5 - X1)(X5 - X2)(X5 - X3))
C...       + .................................... U(X4)
C...         (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5)
C...
C...
C...           + ((X5 - X1)(X5 - X2)(X5 - X3) +
C...
C...              (X5 - X1)(X5 - X2)(X5 - X4) +
C...
C...              (X5 - X1)(X5 - X3)(X5 - X4) +
C...
C...              (X5 - X2)(X5 - X3)(X5 - X4))
C...         .................................... U(X5)
C...         (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4)
C...
C...
C...  EQUATIONS (3) TO (7) ARE THE FIVE-POINT DIFFERENTIATION FORMULAS
C...  FOR A VARIABLE GRID.
C...
C...  NOTE THAT THE FACTORS WHICH MULTIPLY U(X1), U(X2), U(X3), U(X4)
C...  AND U(X5) IN EQUATIONS (3) TO (7) ARE JUST CONSTANTS ONCE THE
C...  VALUES OF X1, X2, X3, X4 AND X5 ARE DEFINED.  THUS, EQUATIONS
C...  (3) TO (7) CAN BE WRITTEN IN THE GENERAL FORM
C...
C...     UX(XI) = C1*U(X1) + C2*U(X2) + C3*U(X3) + C4*U(X4) + C5*U(X5)
C...
C...  WHERE XI = X1, X2, X3, X4, X5 FOR EQUATIONS (3) TO (7) RESPECT-
C...  IVELY, AND C1, C2, C3, C4 AND C5 ARE THE WEIGHTING COEFFICIENTS
C...  CALCULATED IN DSS32A AND PASSED TO DSS032 AS THE THIRD ARGUMENT
C...  (ARRAY C).
C...
C...  IF THE FOLLOWING SUBSTUTIONS FOR A UNIFORM GRID ARE MADE IN EQUA-
C...  TIONS (3) TO (7)
C...
C...
C...     DX = (X2 - X1) = (X3 - X2) = (X4 - X3) = (X5 - X4)
C...
C...    2DX = (X3 - X1) = (X4 - X2) = (X5 - X3)
C...
C...    3DX = (X4 - X1) = (X5 - X2)
C...
C...    4DX = (X5 - X1)
C...
C...    -DX = (X1 - X2) = (X2 - X3) = (X3 - X4) = (X4 - X5)
C...
C...   -2DX = (X1 - X3) = (X2 - X4) = (X3 - X5)
C...
C...   -3DX = (X1 - X4) = (X2 - X5)
C...
C...   -4DX = (X1 - X5)
C...
C...
C...
C...  EQUATIONS (3) TO (7) REDUCE TO THE FOLLOWING FIVE-POINT DIFFER-
C...  TIATION FORMULAE FOR A UNIFORM GRID
C...
C...
C...      UX(X1) = ((-1)(-2)(-3) +
C...
C...                (-1)(-2)(-4) +
C...
C...                (-1)(-3)(-4) +
C...
C...                (-2)(-3)(-4))
C...              ................. U(X1)/DX
C...               (-1)(-2)(-3)(-4)
C...
C...
C...                ((-2)(-3)(-4))
C...            + ................. U(X2)/DX
C...               ( 1)(-1)(-2)(-3)
C...
C...
C...                ((-1)(-3)(-4))
C...            + ................. U(X3)/DX
C...               ( 2)( 1)(-1)(-2)
C...
C...
C...                ((-1)(-2)(-4))
C...            + ................. U(X4)/DX
C...               ( 3)( 2)( 1)(-1)
C...
C...
C...                ((-1)(-2)(-3))
C...            + ................. U(X5)/DX
C...               ( 4)( 3)( 2)( 1)
C...
C...
C...  OR NUMERICALLY EVALUATING THE COEFFICIENTS OF THE U(X1) TO U(X5)
C...  TERMS (WITH 4F = 1*2*3*4 = 4 FACTORIAL)
C...
C...
C...      UX(X1) = (-50/4F)U(X1)/DX
C...
C...             + ( 96/4F)U(X2)/DX
C...
C...             + (-72/4F)U(X3)/DX                                  (8)
C...
C...             + ( 32/4F)U(X4)/DX
C...
C...             + ( -6/4F)U(X5)/DX
C...
C...
C...  SIMILARLY, FROM EQUATIONS (4) TO (7), THE FOLLOWING DIFFERENTIA-
C...  TION FORMULAS FOR A UNIFORM GRID RESULT
C...
C...
C...      UX(X2) =  ((-1)(-2)(-3))
C...              ................. U(X1)/DX
C...               (-1)(-2)(-3)(-4)
C...
C...
C...             + (( 1)(-1)(-2) +
C...
C...                ( 1)(-1)(-3) +
C...
C...                ( 1)(-2)(-3) +
C...
C...                (-1)(-2)(-3))
C...              ................. U(X2)/DX
C...               ( 1)(-1)(-2)(-3)
C...
C...
C...                (( 1)(-2)(-3))
C...            + ................. U(X3)/DX
C...               ( 2)( 1)(-1)(-2)
C...
C...
C...                (( 1)(-1)(-3))
C...            + ................. U(X4)/DX
C...               ( 3)( 2)( 1)(-1)
C...
C...
C...                (( 1)(-1)(-2))
C...            + ................. U(X5)/DX
C...               ( 4)( 3)( 2)( 1)
C...
C...
C...      UX(X2) = ( -6/4F)U(X1)/DX
C...
C...             + (-20/4F)U(X2)/DX
C...
C...             + ( 36/4F)U(X3)/DX                                  (9)
C...
C...             + (-12/4F)U(X4)/DX
C...
C...             + (  2/4F)U(X5)/DX
C...
C...
C...
C...      UX(X3) =  (( 1)(-1)(-2))
C...            + ................. U(X1)/DX
C...               (-1)(-2)(-3)(-4)
C...
C...
C...                (( 2)(-1)(-2))
C...            + ................. U(X2)/DX
C...               ( 1)(-1)(-2)(-3)
C...
C...
C...             + (( 2)( 1)(-1) +
C...
C...                ( 2)( 1)(-2) +
C...
C...                ( 2)(-1)(-2) +
C...
C...                ( 1)(-1)(-2))
C...              ................. U(X3)/DX
C...               ( 2)( 1)(-1)(-2)
C...
C...
C...                (( 2)( 1)(-2))
C...            + ................. U(X4)/DX
C...               ( 3)( 2)( 1)(-1)
C...
C...
C...                (( 2)( 1)(-1))
C...            + ................. U(X5)/DX
C...               ( 4)( 3)( 2)( 1)
C...
C...
C...      UX(X3) = (  2/4F)U(X1)/DX
C...
C...             + (-16/4F)U(X2)/DX
C...
C...             + (  0/4F)U(X3)/DX                                 (10)
C...
C...             + ( 16/4F)U(X4)/DX
C...
C...             + ( -2/4F)U(X5)/DX
C...
C...
C...
C...      UX(X4) =  (( 2)( 1)(-1))
C...              ................. U(X1)/DX
C...               (-1)(-2)(-3)(-4)
C...
C...
C...                (( 3)( 1)(-1))
C...            + ................. U(X2)/DX
C...               ( 1)(-1)(-2)(-3)
C...
C...
C...                (( 3)( 2)(-1))
C...            + ................. U(X3)/DX
C...               ( 2)( 1)(-1)(-2)
C...
C...
C...             + (( 3)( 2)( 1) +
C...
C...                ( 3)( 2)(-1) +
C...
C...                ( 3)( 1)(-1) +
C...
C...                ( 2)( 1)(-1))
C...              ................. U(X4)/DX
C...               ( 3)( 2)( 1)(-1)
C...
C...
C...                (( 3)( 2)( 1))
C...            + ................. U(X5)/DX
C...               ( 4)( 3)( 2)( 1)
C...
C...
C...      UX(X4) = ( -2/4F)U(X1)/DX
C...
C...             + ( 12/4F)U(X2)/DX
C...
C...             + (-36/4F)U(X3)/DX                                 (11)
C...
C...             + ( 20/4F)U(X4)/DX
C...
C...             + (  6/4F)U(X5)/DX
C...
C...
C...
C...      UX(X5) =  (( 3)( 2)( 1))
C...              ................. U(X1)/DX
C...               (-1)(-2)(-3)(-4)
C...
C...
C...                (( 4)( 2)( 1))
C...            + ................. U(X2)/DX
C...               ( 1)(-1)(-2)(-3)
C...
C...
C...                (( 4)( 3)( 1))
C...            + ................. U(X3)/DX
C...               ( 2)( 1)(-1)(-2)
C...
C...
C...                (( 4)( 3)( 2))
C...            + ................. U(X4)/DX
C...               ( 3)( 2)( 1)(-1)
C...
C...
C...             + (( 4)( 3)( 2) +
C...
C...                ( 4)( 3)( 1) +
C...
C...                ( 4)( 2)( 1) +
C...
C...                ( 3)( 2)( 1))
C...              ................. U(X5)/DX
C...               ( 4)( 3)( 2)( 1)
C...
C...
C...      UX(X5) = (  6/4F)U(X1)/DX
C...
C...             + (-32/4F)U(X2)/DX
C...
C...             + ( 72/4F)U(X3)/DX                                 (12)
C...
C...             + (-96/4F)U(X4)/DX
C...
C...             + ( 50/4F)U(X5)/DX
C...
C...
C...  EQUATIONS (8) TO (12) CAN BE SUMMARIZED IN MATRIX FORM AS
C...
C...
C...     ..    ..           ..                           .. ..   ..
C...     .UX(X1).           . -50    96   -72    32    -6 . .U(X1).
C...     .      .           .                             . .     .
C...     .UX(X2).           .  -6   -20    36   -12     2 . .U(X2).
C...     .      .           .                             . .     .
C...     .UX(X3). = (1/4FDX).   2   -16     0    16    -2 .*.U(X3). (13)
C...     .      .           .                             . .     .
C...     .UX(X4).           .  -2    12   -36    20     6 . .U(X4).
C...     .      .           .                             . .     .
C...     .UX(X5).           .   6   -32    72   -96    50 . .U(X5).
C...     ..    ..           ..                           .. ..   ..
C...
C...
C...  THE 5 X 5 COEFFICIENT MATRIX, WITH THE FACTOR (1/4F), IS A BICKLEY
C...  MATRIX FOR N = 4, M = 1, P = 0, 1, 2, 3, 4 (BICKLEY, W. G., FORMU-
C...  LAE FOR NUMERICAL DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION      C,   SUM1,   SUM2,   SUM3,   SUM4,   SUM5,
     1                   TERM,  TERM1,  TERM2,  TERM3,  TERM4,  TERM5,
     2                      X
C...
C...  ABSOLUTE DIMENSION THE ARRAYS USED BY SUBROUTINE DSS32A
      DIMENSION X(5), C(5)
C...
C...  SELECT ONE OF THE DIFFERENTATION FORMULAS, EQUATIONS (3) TO (7)
      GO TO(1,2,3,4,5),I
C...
C...  BEGIN THE CALCULATION OF THE WEIGHTING COEFFICIENTS FOR FIVE-
C...  POINT NUMERICAL DIFFERENTIATION OVER A VARIABLE GRID
C...
C...  ******************************************************************
C...
C...  EQUATION (3)
C...
C...  COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3),
C...  U(X4), U(X5) TERMS IN EQUATION (3) BY DO LOOP 18, STARTING WITH
C...  THE DENOMINATOR OF EACH COEFFICIENT
1     DO 18 I=1,5
      TERM=1.D+00
      DO 11 J=1,5
      IF(J.EQ.I)GO TO 11
      TERM=TERM/(X(I)-X(J))
11    CONTINUE
      IF(I.EQ.1)GO TO 13
C...
C...  COMPLETE THE COEFFICIENTS OF THE U(X2), U(X3), U(X4), U(X5) TERMS
C...  IN EQUATION (3) BY DO LOOP 12, ENDING WITH THE NUMERATOR OF EACH
C...  COEFFICIENT
      DO 12 J=1,5
      IF(J.EQ.1)GO TO 12
      IF(J.EQ.I)GO TO 12
      TERM=(X(1)-X(J))*TERM
12    CONTINUE
      C(I)=TERM
      GO TO 18
C...
C...  COMPLETE THE COEFFICIENT OF THE U(X1) TERM IN EQUATION (3) BY DO
C...  LOOP 15
13    SUM1=0.D+00
      DO 15 JSUM=1,5
      IF(JSUM.EQ.1)GO TO 15
      TERM1=1.D+00
      DO 14 J=1,5
      IF(J.EQ.1)GO TO 14
      IF(J.EQ.JSUM)GO TO 14
      TERM1=TERM1*(X(1)-X(J))
14    CONTINUE
      SUM1=SUM1+TERM1
15    CONTINUE
      C(I)=SUM1*TERM
18    CONTINUE
      RETURN
C...
C...  ******************************************************************
C...
C...  EQUATION (4)
C...
C...  COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3),
C...  U(X4), U(X5) TERMS IN EQUATION (4) BY DO LOOP 28, STARTING WITH
C...  THE DENOMINATOR OF EACH COEFFICIENT
2     DO 28 I=1,5
      TERM=1.D+00
      DO 21 J=1,5
      IF(J.EQ.I)GO TO 21
      TERM=TERM/(X(I)-X(J))
21    CONTINUE
      IF(I.EQ.2)GO TO 23
C...
C...  COMPLETE THE COEFFICIENTS OF THE U(X1), U(X3), U(X4), U(X5) TERMS
C...  IN EQUATION (4) BY DO LOOP 22, ENDING WITH THE NUMERATOR OF EACH
C...  COEFFICIENT
      DO 22 J=1,5
      IF(J.EQ.2)GO TO 22
      IF(J.EQ.I)GO TO 22
      TERM=(X(2)-X(J))*TERM
22    CONTINUE
      C(I)=TERM
      GO TO 28
C...
C...  COMPLETE THE COEFFICIENT OF THE U(X2) TERM IN EQUATION (4) BY DO
C...  LOOP 25
23    SUM2=0.D+00
      DO 25 JSUM=1,5
      IF(JSUM.EQ.2)GO TO 25
      TERM2=1.D+00
      DO 24 J=1,5
      IF(J.EQ.2)GO TO 24
      IF(J.EQ.JSUM)GO TO 24
      TERM2=TERM2*(X(2)-X(J))
24    CONTINUE
      SUM2=SUM2+TERM2
25    CONTINUE
      C(I)=SUM2*TERM
28    CONTINUE
      RETURN
C...
C...  ******************************************************************
C...
C...  EQUATION (5)
C...
C...  COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3),
C...  U(X4), U(X5) TERMS IN EQUATION (5) BY DO LOOP 38, STARTING WITH
C...  THE DENOMINATOR OF EACH COEFFICIENT
3     DO 38 I=1,5
      TERM=1.D+00
      DO 31 J=1,5
      IF(J.EQ.I)GO TO 31
      TERM=TERM/(X(I)-X(J))
31    CONTINUE
      IF(I.EQ.3)GO TO 33
C...
C...  COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X4), U(X5) TERMS
C...  IN EQUATION (5) BY DO LOOP 32, ENDING WITH THE NUMERATOR OF EACH
C...  COEFFICIENT
      DO 32 J=1,5
      IF(J.EQ.3)GO TO 32
      IF(J.EQ.I)GO TO 32
      TERM=(X(3)-X(J))*TERM
32    CONTINUE
      C(I)=TERM
      GO TO 38
C...
C...  COMPLETE THE COEFFICIENT OF THE U(X3) TERM IN EQUATION (5) BY DO
C...  LOOP 35
33    SUM3=0.D+00
      DO 35 JSUM=1,5
      IF(JSUM.EQ.3)GO TO 35
      TERM3=1.D+00
      DO 34 J=1,5
      IF(J.EQ.3)GO TO 34
      IF(J.EQ.JSUM)GO TO 34
      TERM3=TERM3*(X(3)-X(J))
34    CONTINUE
      SUM3=SUM3+TERM3
35    CONTINUE
      C(I)=SUM3*TERM
38    CONTINUE
      RETURN
C...
C...  ******************************************************************
C...
C...  EQUATION (6)
C...
C...  COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3),
C...  U(X4), U(X5) TERMS IN EQUATION (6) BY DO LOOP 48, STARTING WITH
C...  THE DENOMINATOR OF EACH COEFFICIENT
4     DO 48 I=1,5
      TERM=1.D+00
      DO 41 J=1,5
      IF(J.EQ.I)GO TO 41
      TERM=TERM/(X(I)-X(J))
41    CONTINUE
      IF(I.EQ.4)GO TO 43
C...
C...  COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X3), U(X5) TERMS
C...  IN EQUATION (6) BY DO LOOP 42, ENDING WITH THE NUMERATOR OF EACH
C...  COEFFICIENT
      DO 42 J=1,5
      IF(J.EQ.4)GO TO 42
      IF(J.EQ.I)GO TO 42
      TERM=(X(4)-X(J))*TERM
42    CONTINUE
      C(I)=TERM
      GO TO 48
C...
C...  COMPLETE THE COEFFICIENT OF THE U(X4) TERM IN EQUATION (6) BY DO
C...  LOOP 45
43    SUM4=0.D+00
      DO 45 JSUM=1,5
      IF(JSUM.EQ.4)GO TO 45
      TERM4=1.D+00
      DO 44 J=1,5
      IF(J.EQ.4)GO TO 44
      IF(J.EQ.JSUM)GO TO 44
      TERM4=TERM4*(X(4)-X(J))
44    CONTINUE
      SUM4=SUM4+TERM4
45    CONTINUE
      C(I)=SUM4*TERM
48    CONTINUE
      RETURN
C...
C...  ******************************************************************
C...
C...  EQUATION (7)
C...
C...  COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3),
C...  U(X4), U(X5) TERMS IN EQUATION (7) BY DO LOOP 58, STARTING WITH
C...  THE DENOMINATOR OF EACH COEFFICIENT
5     DO 58 I=1,5
      TERM=1.D+00
      DO 51 J=1,5
      IF(J.EQ.I)GO TO 51
      TERM=TERM/(X(I)-X(J))
51    CONTINUE
      IF(I.EQ.5)GO TO 53
C...
C...  COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X3), U(X4) TERMS
C...  IN EQUATION (7) BY DO LOOP 52, ENDING ITH THE NUMERATOR OF EACH
C...  COEFFICIENT
      DO 52 J=1,5
      IF(J.EQ.5)GO TO 52
      IF(J.EQ.I)GO TO 52
      TERM=(X(5)-X(J))*TERM
52    CONTINUE
      C(I)=TERM
      GO TO 58
C...
C...  COMPLETE THE COEFFICIENT OF THE U(X5) TERM IN EQUATION (7) BY DO
C...  LOOP 55
53    SUM5=0.D+00
      DO 55 JSUM=1,5
      IF(JSUM.EQ.5)GO TO 55
      TERM5=1.D+00
      DO 54 J=1,5
      IF(J.EQ.5)GO TO 54
      IF(J.EQ.JSUM)GO TO 54
      TERM5=TERM5*(X(5)-X(J))
54    CONTINUE
      SUM5=SUM5+TERM5
55    CONTINUE
      C(I)=SUM5*TERM
58    CONTINUE
      RETURN
C...
C...  ******************************************************************
C...
      END
*DECK DSS034
      SUBROUTINE DSS034(XL,XU,N1,N2,ND,U2D,UX2D,V)
C...
C...  SUBROUTINE DSS034 COMPUTES A PARTIAL DERIVATIVE OVER A TWO-
C...  DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE-
C...  POINT BIASED UPWIND APPROXIMATIONS.  IT IS INTENDED PRIMARILY
C...  FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION
C...  OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN TWO DIMENSIONS.
C...  DSS034 IS RECOMMENDED PARTICULARLY FOR CONVECTIVE-DIFFUSION
C...  PROBLEMS WHICH REQUIRE COMBINATIONS OF BIASED UPWIND AND CENTERED
C...  APPROXIMATIONS FOR FIRST AND SECOND-ORDER SPATIAL DERIVATIVES IN
C...  PDES.
C...
C...  ARGUMENT LIST
C...
C...     XL        LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     XU        UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     N1        NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     N2        NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     ND        NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE
C...               PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     U2D       TWO-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT VARI-
C...               ABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (INPUT)
C...
C...     UX2D      TWO-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI-
C...               VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (OUTPUT)
C...
C...     V         VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED
C...               OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE
C...               PARTIAL DERIVATIVE.  V EQ 0 CALLS THE FIVE-POINT
C...               CENTERED APPROXIMATION.  V NE 0 CALLS THE FIVE-POINT
C...               BIASED UPWIND APPROXIMATION (INPUT)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION   UX1D,   UX2D,    U1D,    U2D,      V,     XL,
     1                     XU
C...
C...  THE FOLLOWING TWO-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U2D) AND ITS PARTIAL DERIVATIVE (UX2D)
      DIMENSION   U2D(N1,N2), UX2D(N1,N2)
C...
C...  THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D).  IN EACH
C...  CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE
C...  OTHER INDEPENDENT VARIABLE VARIES OVER ITS TOTAL INTERVAL.
C...  THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE
C...  ONE-DIMENSIONAL ROUTINES DSS004 AND DSS020.
C...
C...  NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE-
C...  FORE HAVE TO BE INCREASED IN SIZE.  HOWEVER, WITH A SIZE
C...  OF 51, THE TWO-DIMENSIONAL PROBLEM COULD HAVE A GRID OF
C...  51 X 51 POINTS, THEREBY GENERATING AN APPROXIMATING ODE
C...  SYSTEM WITH A MULTIPLE OF 51 X 51 EQUATIONS, DEPENDING ON
C...  THE NUMBER OF SIMULTANEOUS PDES.  THIS IS A VERY LARGE ODE
C...  PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING
C...  IS CONSIDERED ADEQUATE FOR MOST PROBLEMS.
      DIMENSION      U1D(51),    UX1D(51)
C...
C...  GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED
C...  WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE
      IF(ND.EQ.2)GO TO 2
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N1 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 10, 11 AND 12
      DO 10 J=1,N2
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 11 I=1,N1
      U1D(I)=U2D(I,J)
11    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS020(XL,XU,N1,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UX2D
      DO 12 I=1,N1
      UX2D(I,J)=UX1D(I)
12    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE SECOND INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE SECOND INDEPENDENT VARIABLE
10    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D
      RETURN
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N2 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 GRID POINTS VIA NESTED DO LOOPS 20 AND 21
2     DO 20 I=1,N1
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 21 J=1,N2
      U1D(J)=U2D(I,J)
21    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS020(XL,XU,N2,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE TWO-DIMENSIONAL ARRAY UX2D
      DO 22 J=1,N2
      UX2D(I,J)=UX1D(J)
22    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE FIRST INDE-
C...  PENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE CALCULATION FOR
C...  THE NEXT VALUE OF THE FIRST INDEPENDENT VARIABLE
20    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL
C...  DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D
      RETURN
      END
*DECK DSS036
      SUBROUTINE DSS036(XL,XU,N1,N2,N3,ND,U3D,UX3D,V)
C...
C...  SUBROUTINE DSS036 COMPUTES A PARTIAL DERIVATIVE OVER A THREE-
C...  DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE-
C...  POINT BIASED UPWIND APPROXIMATIONS.  IT IS INTENDED PRIMARILY
C...  FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION
C...  OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN THREE DIMENSIONS.
C...
C...  SUBROUTINE DSS036 IS CALLED IN ESSENTIALLY THE SAME WAY AS
C...  SUBROUTINE DSS034.  THE ONLY DIFFERENCE IS AN ADDITIONAL ARGU-
C...  MENT, N3, TO DEFINE THE NUMBER OF GRID POINTS IN THE THIRD
C...  DIMENSION.  THE COMMENTS IN DSS034 SHOULD THEREFORE BE USEFUL
C...  IN UNDERSTANDING THE OPERATION OF DSS036.  IN PARTICULAR,
C...  DSS036 CALLS SUBROUTINES DSS004 AND DSS020 TO IMPLEMENT THE
C...  FIVE-POINT CENTERED APPROXIMATION AND FIVE-POINT BIASED UPWIND
C...  APPROXIMATION OF THE PARTIAL DERIAVTIVE, RESPECTIVELY.
C...
C...  ARGUMENT LIST
C...
C...     XL        LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     XU        UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH
C...               THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     N1        NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     N2        NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     N3        NUMBER OF GRID POINTS FOR THE THIRD INDEPENDENT
C...               VARIABLE (INPUT)
C...
C...     ND        NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE
C...               PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT)
C...
C...     U3D       THREE-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT
C...               VARIABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT
C...               TO INDEPENDENT VARIABLE ND (INPUT)
C...
C...     UX3D      THREE-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI-
C...               VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO
C...               INDEPENDENT VARIABLE ND (OUTPUT)
C...
C...     V         VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED
C...               OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE
C...               PARTIAL DERIVATIVE.  V EQ 0 CALLS THE FIVE-POINT
C...               CENTERED APPROXIMATION.  V NE 0 CALLS THE FIVE-POINT
C...               BIASED UPWIND APPROXIMATION (INPUT)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION   UX1D,   UX3D,    U1D,    U3D,      V,     XL,
     1                     XU
C...
C...  THE FOLLOWING THREE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U3D) AND ITS PARTIAL DERIVATIVE (UX3D)
      DIMENSION   U3D(N1,N2,N3), UX3D(N1,N2,N3)
C...
C...  THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT
C...  VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D).  IN EACH
C...  CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE
C...  OTHER TWO INDEPENDENT VARIABLES VARY OVER THEIR TOTAL INTERVALS.
C...  THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE
C...  ONE-DIMENSIONAL ROUTINES DSS004 AND DSS020.
C...
C...  NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE-
C...  FORE HAVE TO BE INCREASED IN SIZE.  HOWEVER, WITH A SIZE
C...  OF 51, THE THREE-DIMENSIONAL PROBLEM COULD HAVE A GRID OF
C...  51 X 51 X 51 POINTS, THEREBY GENERATING AN APPROXIMATING ODE
C...  SYSTEM WITH A MULTIPLE OF 51 X 51 X 51 EQUATIONS, DEPENDING ON
C...  THE NUMBER OF SIMULTANEOUS PDES.  THIS IS A VERY LARGE ODE
C...  PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING
C...  IS CONSIDERED ADEQUATE FOR MOST PROBLEMS.
      DIMENSION      U1D(51),    UX1D(51)
C...
C...  GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED
C...  WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE
      IF(ND.EQ.2)GO TO 2
C...
C...  GO TO STATEMENT 3 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED
C...  WITH RESPECT TO THE THIRD INDEPENDENT VARIABLE
      IF(ND.EQ.3)GO TO 3
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N1 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 X N3 GRID POINTS VIA NESTED DO LOOPS 10, 11, 12 AND 13
      DO 10 J=1,N2
      DO 11 K=1,N3
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 12 I=1,N1
      U1D(I)=U3D(I,J,K)
12    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS020(XL,XU,N1,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE THREE-DIMENSIONAL ARRAY UX3D
      DO 13 I=1,N1
      UX3D(I,J,K)=UX1D(I)
13    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE SECOND AND
C...  THIRD INDEPENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE
C...  CALCULATION FOR THE OTHER VALUES OF THE SECOND AND THIRD
C...  INDEPENDENT VARIABLES
11    CONTINUE
10    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 X N3 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE
C...  PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D
      RETURN
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N2 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 X N3 GRID POINTS VIA NESTED DO LOOPS 20, 21, 22 AND 23
2     DO 20 I=1,N1
      DO 21 K=1,N3
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 22 J=1,N2
      U1D(J)=U3D(I,J,K)
22    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS020(XL,XU,N2,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE THREE-DIMENSIONAL ARRAY UX3D
      DO 23 J=1,N2
      UX3D(I,J,K)=UX1D(J)
23    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE FIRST AND
C...  THIRD INDEPENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE
C...  CALCULATION FOR THE OTHER VALUES OF THE FIRST AND THIRD
C...  INDEPENDENT VARIABLES
21    CONTINUE
20    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 X N3 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE
C...  PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D
      RETURN
C...
C...  ******************************************************************
C...
C...  THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE
C...  THIRD INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING
C...  OF N3 GRID POINTS.  COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X
C...  N2 X N3 GRID POINTS VIA NESTED DO LOOPS 30, 31, 32 AND 33
3     DO 30 I=1,N1
      DO 31 J=1,N2
C...
C...  TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D
C...  TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND
C...  DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE
      DO 32 K=1,N3
      U1D(K)=U3D(I,J,K)
32    CONTINUE
C...
C...  IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE
C...  PARTIAL DERIVATIVE
      IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N3,U1D,UX1D)
C...
C...  IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR
C...  THE PARTIAL DERIVATIVE
      IF(V.NE.0.D+00)CALL DSS020(XL,XU,N3,U1D,UX1D,V)
C...
C...  RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D
C...  TO THE THREE-DIMENSIONAL ARRAY UX3D
      DO 33 K=1,N3
      UX3D(I,J,K)=UX1D(K)
33    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE FIRST AND
C...  SECOND INDEPENDENT VARIABLE HAS BEEN CALCULATED.  REPEAT THE
C...  CALCULATION FOR THE OTHER VALUES OF THE FIRST AND SECOND
C...  INDEPENDENT VARIABLES
31    CONTINUE
30    CONTINUE
C...
C...  THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X
C...  N2 X N3 GRID.  THEREFORE RETURN TO THE CALLING PROGRAM WITH THE
C...  PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D
      RETURN
      END
*DECK DSS038
      SUBROUTINE DSS038(N,X,U,UX,B,C,D)
C...
C...  SUBROUTINE DSS038 IS A NUMERICAL DIFFERENTAITION ROUTINE BASED ON
C...  THE CUBIC SPLINE OF FORSYTHE, G. E., M. A. MALCOLM, AND C. B.
C...  MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-
C...  HALL, ENGLEWOOD CLIFFS, N.J., 1977, PP. 76-79.
C...
C...  ARGUMENT LIST
C...
C...     N      NUMBER OF SPATIAL GRID POINTS (INPUT)
C...
C...     X      ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE
C...            SPATIAL (BOUNDARY-VALUE) INDEPENDENT VARIABLE (INPUT)
C...
C...     U      ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE
C...            DEPENDENT VARIABLE (INPUT)
C...
C...     UX     ONE-DIMENSIONAL ARRAY CONTAINING THE DERIVATIVE OF U
C...            WITH RESPECT TO X (OUTPUT)
C...
C...     B,C,D  ONE-DIMENSIONAL ARRAYS CONTAINING THE CUBIC SPLINE CO-
C...            EFFICIENTS (OUTPUT)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION      B,      C,      D,      U,     UX,      X
C...
C...  VARIABLE DIMENSION THE ARRAYS PASSED AS ARGUMENTS
      DIMENSION   X(N),  U(N), UX(N),  B(N),  C(N),  D(N)
C...
C...  CALL SUBROUTINE SPLINE TO COMPUTE THE SPLINE COEFFICIENTS
      CALL SPLINE(N,X,U,B,C,D)
C...
C...  THE NUMERICAL DERIVATIVE IS OBTAINED BY DIFFERENTIATION OF THE
C...  CUBIC SPLINE
C...
C...     S(X) = U(I) + B(I)*(X - X(I)) + C(I)*(X - X(I))**2
C...
C...                 + D(I)*(X - X(I))**3
C...  THUS
C...
C...     UX(I) = DS(X(I))/DX = B(I)
C...
      DO 1 I=1,N
      UX(I)=B(I)
1     CONTINUE
      RETURN
      END
*DECK DSS040
      SUBROUTINE DSS040(N,X,U,UX,UXX,B,C,D)
C...
C...  SUBROUTINE DSS040 IS A NUMERICAL DIFFERENTAITION ROUTINE BASED ON
C...  THE CUBIC SPLINE OF FORSYTHE, G. E., M. A. MALCOLM, AND C. B.
C...  MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-
C...  HALL, ENGLEWOOD CLIFFS, N.J., 1977, PP. 76-79.
C...
C...  ARGUMENT LIST
C...
C...     N      NUMBER OF SPATIAL GRID POINTS (INPUT)
C...
C...     X      ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE
C...            SPATIAL (BOUNDARY-VALUE) INDEPENDENT VARIABLE (INPUT)
C...
C...     U      ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE
C...            DEPENDENT VARIABLE (INPUT)
C...
C...     UX     ONE-DIMENSIONAL ARRAY CONTAINING THE FIRST DERIVATIVE
C...            OF U WITH RESPECT TO X (OUTPUT)
C...
C...     UXX    ONE-DIMENSIONAL ARRAY CONTAINING THE SECOND DERIVATIVE
C...            OF U WITH RESPECT TO X (OUTPUT)
C...
C...     B,C,D  ONE-DIMENSIONAL ARRAYS CONTAINING THE CUBIC SPLINE CO-
C...            EFFICIENTS (OUTPUT)
C...
C...  TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION      B,      C,      D,      U,     UX,    UXX,
     1                      X
C...
C...  VARIABLE DIMENSION THE ARRAYS PASSED AS ARGUMENTS
      DIMENSION   X(N),  U(N), UX(N),UXX(N),  B(N),  C(N),  D(N)
C...
C...  CALL SUBROUTINE SPLINE TO COMPUTE THE SPLINE COEFFICIENTS
      CALL SPLINE(N,X,U,B,C,D)
C...
C...  THE NUMERICAL DERIVATIVES, UX AND UXX, ARE COMPUTED BY DIFFERENT-
C...  IATION OF THE CUBIC SPLINE
C...
C...     S(X) = U(I) + B(I)*(X - X(I)) + C(I)*(X - X(I))**2
C...
C...                 + D(I)*(X - X(I))**3
C...  THUS
C...
C...     UX(I) = DS(X(I))/DX = B(I)
C...
C...     UXX(I) = D(DS(X(I))/DX)/DX = 2*C(I)
C...
      DO 1 I=1,N
      UX(I)=B(I)
      UXX(I)=2.D+00*C(I)
1     CONTINUE
      RETURN
      END
*DECK SPLINE
      SUBROUTINE SPLINE (N, X, Y, B, C, D)
C...
C...  THIS SUBROUTINE APPEARS IN
C...
C...     FORSYTHE, G. E., M. A. MALCOLM AND C. B. MOLER, COMPUTER
C...     METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-HALL, INC.,
C...     ENGLEWOOD CLIFFS, N.J., 1977, 76-78
C...
      INTEGER N
      DOUBLE PRECISION X(N), Y(N), B(N), C(N), D(N)
C
C  THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED
C  FOR A CUBIC INTERPOLATING SPLINE
C
C    S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3
C
C    FOR  X(I) .LE. X .LE. X(I+1)
C
C  INPUT..
C
C    N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2)
C    X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER
C    Y = THE ORDINATES OF THE KNOTS
C
C  OUTPUT..
C
C    B, C, D  = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE.
C
C  USING  P  TO DENOTE DIFFERENTIATION,
C
C    Y(I) = S(X(I))
C    B(I) = SP(X(I))
C    C(I) = SPP(X(I))/2
C    D(I) = SPPP(X(I))/6  (DERIVATIVE FROM THE RIGHT)
C
C  THE ACCOMPANYING FUNCTION SUBPROGRAM  SEVAL  CAN BE USED
C  TO EVALUATE THE SPLINE.
C
C
      INTEGER NM1, IB, I
      DOUBLE PRECISION T
C
      NM1 = N-1
      IF ( N .LT. 2 ) RETURN
      IF ( N .LT. 3 ) GO TO 50
C
C  SET UP TRIDIAGONAL SYSTEM
C
C  B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE.
C
      D(1) = X(2) - X(1)
      C(2) = (Y(2) - Y(1))/D(1)
      DO 10 I = 2, NM1
         D(I) = X(I+1) - X(I)
         B(I) = 2.D+00*(D(I-1) + D(I))
         C(I+1) = (Y(I+1) - Y(I))/D(I)
         C(I) = C(I+1) - C(I)
   10 CONTINUE
C
C  END CONDITIONS.  THIRD DERIVATIVES AT  X(1)  AND  X(N)
C  OBTAINED FROM DIVIDED DIFFERENCES
C
      B(1) = -D(1)
      B(N) = -D(N-1)
      C(1) = 0.D+00
      C(N) = 0.D+00
      IF ( N .EQ. 3 ) GO TO 15
      C(1) = C(3)/(X(4)-X(2)) - C(2)/(X(3)-X(1))
      C(N) = C(N-1)/(X(N)-X(N-2)) - C(N-2)/(X(N-1)-X(N-3))
      C(1) = C(1)*D(1)**2/(X(4)-X(1))
      C(N) = -C(N)*D(N-1)**2/(X(N)-X(N-3))
C
C  FORWARD ELIMINATION
C
   15 DO 20 I = 2, N
         T = D(I-1)/B(I-1)
         B(I) = B(I) - T*D(I-1)
         C(I) = C(I) - T*C(I-1)
   20 CONTINUE
C
C  BACK SUBSTITUTION
C
      C(N) = C(N)/B(N)
      DO 30 IB = 1, NM1
         I = N-IB
         C(I) = (C(I) - D(I)*C(I+1))/B(I)
   30 CONTINUE
C
C  C(I) IS NOW THE SIGMA(I) OF THE TEXT
C
C  COMPUTE POLYNOMIAL COEFFICIENTS
C
      B(N) = (Y(N) - Y(NM1))/D(NM1) + D(NM1)*(C(NM1) + 2.D+00*C(N))
      DO 40 I = 1, NM1
         B(I) = (Y(I+1) - Y(I))/D(I) - D(I)*(C(I+1) + 2.D+00*C(I))
         D(I) = (C(I+1) - C(I))/D(I)
         C(I) = 3.D+00*C(I)
   40 CONTINUE
      C(N) = 3.D+00*C(N)
      D(N) = D(N-1)
      RETURN
C
   50 B(1) = (Y(2)-Y(1))/(X(2)-X(1))
      C(1) = 0.D+00
      D(1) = 0.D+00
      B(2) = B(1)
      C(2) = 0.D+00
      D(2) = 0.D+00
      RETURN
      END
*DECK DSS042
      SUBROUTINE DSS042(XL,XU,N,U,UX,UXX,NL,NU)
C...
C...  SUBROUTINE DSS042 COMPUTES A SECOND-ORDER APPROXIMATION OF A
C...  SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE
C...  AT THE BOUNDARY.
C...
C...     DOUBLE PRECISION
C...
C...  ARGUMENT LIST
C...
C...     XL      LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
C...
C...     XU      RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
C...
C...     N       NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END
C...             POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE
C...             DIFFERENTIATED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U.
C...             THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN
C...             NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU,
C...             DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE-
C...             SCRIPTION OF NL AND NU BELOW)
C...
C...     UXX     ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U
C...             (OUTPUT)
C...
C...     NL      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
C...             X = XL (INPUT).  THE ALLOWABLE VALUES ARE
C...
C...                1 - DIRICHLET BOUNDARY CONDITION AT X = XL
C...                    (UX(1) IS NOT USED)
C...
C...                2 - NEUMANN BOUNDARY CONDITION AT X = XL
C...                    (UX(1) IS USED)
C...
C...     NU      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
C...             X = XU (INPUT).  THE ALLOWABLE VALUES ARE
C...
C...                1 - DIRICHLET BOUNDARY CONDITION AT X = XU
C...                    (UX(N) IS NOT USED)
C...
C...                2 - NEUMANN BOUNDARY CONDITION AT X = XU
C...                    (UX(N) IS USED)
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
      DIMENSION U(N), UX(N), UXX(N)
C...
C...  THE FOLLOWING DERIVATION IS FOR A SET OF SECOND-ORDER, FOUR-POINT
C...  APPROXIMATIONS FOR A SECOND DERIVATIVE THAT CAN BE USED AT THE
C...  BOUNDARIES OF A SPATIAL DOMAIN.  THESE APPROXIMATIONS HAVE THE
C...  FEATURES
C...
C...     (1)  ONLY INTERIOR AND BOUNDARY POINTS ARE USED (I.E., NO
C...          FICTITIOUS POINTS ARE USED)
C...
C...     (2)  THE NORMAL DERIVATIVE AT THE BOUNDARY IS INCLUDED AS PART
C...          OF THE APPROXIMATION FOR THE SECOND DERIVATIVE
C...
C...     (3)  APPROXIMATIONS FOR THE BOUNDARY CONDITIONS ARE NOT USED.
C...
C...  THE DERIVATION IS BY PROFESSOR GILBERT A. STENGLE, DEPARTMENT OF
C...  MATHEMATICS, LEHIGH UNIVERSITY, BETHLEHEM, PA 18015, AND WAS DONE
C...  ON DECEMBER 7, 1985.
C...
C...  FOR AN APPROXIMATION AT THE LEFT BOUNDARY, INVOLVING THE POINTS
C...  I, I+1, I+2 AND I+3, CONSIDER THE FOLLOWING TAYLOR SERIES EXPAN-
C...  SIONS
C...
C...                  UX(I)( DX)   UXX(I)( DX)**2   UXXX(I)( DX)**3
C...  U(I+1) = U(I) + ---------- + -------------- + --------------- +...
C...                       1             2                 6
C...
C...
C...                  UX(I)(2DX)   UXX(I)(2DX)**2   UXXX(I)(2DX)**3
C...  U(I+2) = U(I) + ---------- + -------------- + --------------- +...
C...                       1             2                 6
C...
C...  IF WE NOW FORM THE FOLLOWING LINEAR COMBINATION, INVOLVING CON-
C...  STANTS A, B, C AND D TO BE DETERMINED, AND USE THE PRECEDING TWO
C...  TAYLOR SERIES,
C...
C...     A*U(I) + B*UX(I) + C*U(I+1) + D*U(I+2)
C...
C...  WE HAVE
C...
C...     A*U(I) + B*UX(I) + C*U(I+1) + D*U(I+2) =
C...
C...     (A + B + C + D)*U(I) +
C...
C...     (B + DX*C + 2*DX*D)*UX(I) +
C...
C...     (C*(DX**2)/2 + D*((2*DX)**2)/2)*UXX(I) +
C...
C...     (C*(DX**3)/6 + D*((2*DX)**3)/6)*UXXX(I) + O(DX**4)
C...
C...  THE THIRD DERIVATIVE, UXXX(I), CAN BE DROPPED BY TAKING
C...
C...     C = -8*D
C...
C...  THE SECOND DERIVATIVE, UXX(I), CAN BE RETAINED BY TAKING
C...
C...     (DX**2)(C/2 + 2*D) = 1
C...
C...  WHICH, WHEN COMBINED WITH THE PRECEDING RESULT GIVES
C...
C...     D = -1/(2*(DX**2))
C...
C...     C = 4/(DX**2)
C...
C...  THE FIRST DERIVATIVE, UX(I), CAN BE DROPPED BY TAKING
C...
C...     B + DX*C + 2*DX*D = 0
C...
C...  OR
C...
C...     B = -DX*C - 2*DX*D = -4/DX - 2*DX*(-1/(2*(DX**2))) = -3/DX
C...
C...  FINALLY, U(I), CAN BE DROPPED BY TAKING
C...
C...     A = - C - D = 8*D - D = -7*D = -7/(2*(DX**2))
C...
C...  IF WE NOW SOLVE FOR THE DERIVATIVE OF INTEREST, UXX(I),
C...
C...     UXX(I) = -7/(2(DX**2))*U(I) - 3/DX*UX(I)
C...
C...              + 8/(DX**2)*U(I+1) - 1/(2*(DX**2))U(I+2) + O(DX**2)
C...
C...       = (1/(2*(DX**2)))*(-U(I+2) + 8*U(I+1) - 7*U(I) - 6*DX*UX(I))
C...
C...         + O(DX**2)
C...
C...  WHICH IS THE FOUR-POINT, SECOND-ORDER APPROXIMATION FOR THE SECOND
C...  DERIVATIVE, UXX(I), INCLUDING THE FIRST DERIVATIVE, UX(I).
C...
C...  FOUR CHECKS OF THIS APPROXIMATION CAN EASILY BE MADE FOR U(I) =
C...  1, U(I) = X, U(I) = X**2 AND U(I) = X**3
C...
C...     UXX(I) = (1/(2*(DX**2)))*(-1 + 8*1 - 7*1 - 6*DX*0) = 0
C...
C...     UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX) + 8*(X + DX)
C...
C...              -7*X - 6*DX*1) = 0
C...
C...     UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX)**2 + 8*(X + DX)**2
C...
C...            - 7*(X**2) - 6*DX*(2*X))
C...
C...             = (-  X**2 -  4*X*DX - 4*DX**2
C...
C...               + 8*X**2 + 16*X*DX + 8*DX**2
C...
C...               - 7*X**2 - 12*X*DX)/(2*(DX**2)) = 2
C...
C...     UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX)**3 + 8*(X + DX)**3
C...
C...            - 7*(X**3) - 6*DX*(3*X**2))
C...
C...            = (1/(2*(DX**2)))*(- X**3 - 6*DX*X**2 - 12*X*DX**2
C...
C...            - 8*DX**3 + 8*X**3 + 24*DX*X**2 + 24*X*DX**2 + 8*DX**3
C...
C...            - 7*X**3 - 18*DX*X**2)
C...
C...            = (1/(2*(DX**2)))*(12*X*DX**2) = 6*X
C...
C...  THE PRECEDING APPROXIMATION FOR UXX(I) CAN BE APPLIED AT THE
C...  LEFT BOUNDARY VALUE OF X BY TAKING I = 1.  AN APPROXIMATION AT
C...  THE RIGHT BOUNDARY IS OBTAINED BY TAKING DX = -DX AND REVERSING
C...  THE SUBSCRIPTS IN THE PRECEDING APPROXIMATION, WITH I = N
C...
C...     UXX(I)
C...
C...       = (1/(2*(DX**2)))*(-U(I-2) + 8*U(I-1) - 7*U(I) + 6*DX*UX(I))
C...
C...         + O(DX**2)
C...
C...  TO OBTAIN APPROXIMATIONS OF THE SECOND DERVIAVTIVE WHICH DO NOT
C...  INVOLVE THE FIRST DERIVATIVE, WE TAKE AS THE LINEAR COMBINATION
C...
C...     A*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) =
C...
C...  WE HAVE
C...
C...     A*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) =
C...
C...     (A + B + C + D)*U(I)+
C...
C...     (DX*B + 2*DX*C + 4*DX*D)*UX(I)+
C...
C...     (B*(DX**2)/2 + C*((2*DX)**2)/2 + D*((3*DX)**2)/2)*UXX(I) +
C...
C...     (B*(DX**3)/6 + C*((2*DX)**3)/6 + D*((3*DX)**3)/6)*UXX(I) +
C...
C...     O(DX**4)
C...
C...  THE THIRD DERIVATIVE, UXXX(I), CAN BE DROPPED BY TAKING
C...
C...     B + 8*C + 27*D = 0
C...
C...  THE SECOND DERIVATIVE, UXX(I), CAN BE RETAINED BY TAKING
C...
C...     (DX**2)*(B/2 + 2*C + (9/2)*D) = 1
C...
C...  THE FIRST DERIVATIVE CAN BE DROPPED BY TAKING
C...
C...     B + 2*C + 3*D = 0
C...
C...  SOLUTION OF THE PRECEDING EQUATIONS FOR C AND D BY ELIMINATION OF
C...  B GIVES
C...
C...     6*C + 24*D = 0
C...
C...     4*C + 18*D = -2/(DX**2)
C...
C...  THEN, ELIMINATING C, GIVES
C...
C...     (18 - 16)*D = -2/(DX**2)
C...
C...  OR
C...
C...     D = -1/(DX**2)
C...
C...     C = (24/6)/(DX**2) = 4/(DX**2)
C...
C...     B = -8/(DX**2) + 3/(DX**2) = -5/(DX**2)
C...
C...  U(I) CAN BE DROPPED BY TAKING
C...
C...     A + B + C + D = 0
C...
C...  OR
C...
C...     A = (5 - 4 + 1)/(DX**2) = 2/(DX**2)
C...
C...  IF WE NOW SOLVE FOR THE DERIVATIVE OF INTEREST, UXX(I),
C...
C...     UXX(I) = (1/DX**2)*(2*U(I) - 5*U(I+1) + 4*U(I+2) - 1*U(I+3))
C...
C...            + O(DX**2)
C...
C...  WHICH IS THE FOUR-POINT, SECOND-ORDER APPROXIMATION FOR THE SECOND
C...  DERIVATIVE, UXX(I), WITHOUT THE FIRST DERIVATIVE, UX(I).
C...
C...  FOUR CHECKS OF THIS APPROXIMATION CAN EASILY BE MADE FOR U(I) =
C...  1, U(I) = X, U(I) = X**2 AND U(I) = X**3
C...
C...     UXX(I) = (1/DX**2)*(2 - 5 + 4 - 1) = 0
C...
C...     UXX(I) = (1/DX**2)*(2*X - 5*(X + DX) + 4*(X + 2*DX)
C...
C...              - 1*(X + 3*DX)) = 0
C...
C...     UXX(I) = (1/DX**2)*(2*X**2 - 5*(X + DX)**2 + 4*(X + 2*DX)**2
C...
C...              - 1*(X + 3*DX)**2) = 2
C...
C...     UXX(I) = (1/DX**2)*(2*X**3 - 5*(X + DX)**3 + 4*(X + 2*DX)**3
C...
C...            - 1*(X + 3*DX)**3)
C...
C...             = (1/DX**2)*(2*X**3 - 5*X**3 - 15*X*DX**2
C...
C...             - 15*DX*X**2 - 5*DX**3 + 4*X**3 + 24*DX*X**2
C...
C...             + 48*X*DX**2 + 32*DX**3 - X**3 - 9*DX*X**2
C...
C...             - 27*X*DX**2 - 27DX**3)
C...
C...             = (1/DX**2)*(6*X*DX**2) = 6*X
C...
C...  THE PRECEDING APPROXIMATION FOR UXX(I) CAN BE APPLIED AT THE
C...  LEFT BOUNDARY VALUE OF X BY TAKING I = 1.  AN APPROXIMATION AT
C...  THE RIGHT BOUNDARY IS OBTAINED BY TAKING DX = -DX AND REVERSING
C...  THE SUBSCRIPTS IN THE PRECEDING APPROXIMATION, WITH I = N
C...
C...     UXX(I) = (1/DX**2)*(2*U(I) - 5*U(I-1) + 4*U(I-2) - 1*U(I-3))
C...
C...            + O(DX**2)
C...
C...  GRID SPACING
      DX=(XU-XL)/DFLOAT(N-1)
C...
C...  CALCULATE UXX AT THE LEFT BOUNDARY, WITHOUT UX
      IF(NL.EQ.1)THEN
      UXX(1)=((     2.D0)*U(  1)
     1       +(    -5.D0)*U(  2)
     2       +(     4.D0)*U(  3)
     3       +(    -1.D0)*U(  4))/(DX**2)
C...
C...  CALCULATE UXX AT THE LEFT BOUNDARY, INCLUDING UX
      ELSE IF(NL.EQ.2)THEN
      UXX(1)=((    -7.D0)*U(  1)
     1       +(     8.D0)*U(  2)
     2       +(    -1.D0)*U(  3))/(2.D0*DX**2)
     3       +(    -6.D0)*UX( 1) /(2.D0*DX)
      END IF
C...
C...  CALCULATE UXX AT THE RIGHT BOUNDARY, WITHOUT UX
      IF(NU.EQ.1)THEN
      UXX(N)=((     2.D0)*U(N  )
     1       +(    -5.D0)*U(N-1)
     2       +(     4.D0)*U(N-2)
     3       +(    -1.D0)*U(N-3))/(DX**2)
C...
C...  CALCULATE UXX AT THE RIGHT BOUNDARY, INCLUDING UX
      ELSE IF(NU.EQ.2)THEN
      UXX(N)=((    -7.D0)*U(N  )
     1       +(     8.D0)*U(N-1)
     2       +(    -1.D0)*U(N-2))/(2.D0*DX**2)
     3       +(     6.D0)*UX(N ) /(2.D0*DX)
      END IF
C...
C...  CALCULATE UXX AT THE INTERIOR GRID POINTS
      DO 1 I=2,N-1
      UXX(I)=(U(I+1)-2.D0*U(I)+U(I-1))/DX**2
1     CONTINUE
      RETURN
      END
*DECK DSS044
      SUBROUTINE DSS044(XL,XU,N,U,UX,UXX,NL,NU)
C...
C...  SUBROUTINE DSS044 COMPUTES A FOURTH-ORDER APPROXIMATION OF A
C...  SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE
C...  AT THE BOUNDARY.
C...
C...     SINGLE PRECISION
C...
C...  ARGUMENT LIST
C...
C...     XL      LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
C...
C...     XU      RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT)
C...
C...     N       NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END
C...             POINTS (INPUT)
C...
C...     U       ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE
C...             DIFFERENTIATED (INPUT)
C...
C...     UX      ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U.
C...             THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN
C...             NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU,
C...             DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE-
C...             SCRIPTION OF NL AND NU BELOW)
C...
C...     UXX     ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U
C...             (OUTPUT)
C...
C...     NL      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
C...             X = XL (INPUT).  THE ALLOWABLE VALUES ARE
C...
C...                1 - DIRICHLET BOUNDARY CONDITION AT X = XL
C...                    (UX(1) IS NOT USED)
C...
C...                2 - NEUMANN BOUNDARY CONDITION AT X = XL
C...                    (UX(1) IS USED)
C...
C...     NU      INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT
C...             X = XU (INPUT).  THE ALLOWABLE VALUES ARE
C...
C...                1 - DIRICHLET BOUNDARY CONDITION AT X = XU
C...                    (UX(N) IS NOT USED)
C...
C...                2 - NEUMANN BOUNDARY CONDITION AT X = XU
C...                    (UX(N) IS USED)
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
      DIMENSION U(N), UX(N), UXX(N)
C...
C...  THE FOLLOWING DERIVATION WAS COMPLETED BY W. E. SCHIESSER, DEPTS
C...  OF CHE AND CSEE, LEHIGH UNIVERSITY, BETHLEHEM, PA 18015, USA, ON
C...  DECEMBER 15, 1986.  ADDITIONAL DETAILS ARE GIVEN IN SUBROUTINE
C...  DSS042.
C...
C...  ******************************************************************
C...
C...  (1)  UXX AT THE INTERIOR POINTS 3, 4,..., N-2
C...
C...  TO DEVELOP A SET OF FOURTH-ORDER CORRECT DIFFERENTIATION FORMULAS
C...  FOR THE SECOND DERIVATIVE UXX, WE CONSIDER FIRST THE INTERIOR
C...  GRID POINTS AT WHICH A SYMMETRIC FORMULA CAN BE USED.
C...
C...  IF WE CONSIDER A FORMULA OF THE FORM
C...
C...     A*U(I-2) + B*U(I-1) + E*U(I) + C*U(I+1) + D*U(I+2)
C...
C...  TAYLOR SERIES EXPANSIONS OF U(I-2), U(I-1), U(I+1) AND U(I+2)
C...  CAN BE SUBSTITUTED INTO THIS FORMULA.  WE THEN CONSIDER THE
C...  LINEAR ALBEGRAIC EQUATIONS RELATING A, B, C AND D WHICH WILL
C...  RETAIN CERTAIN TERMS, I.E., UXX, AND DROP OTHERS, E.G., UXXX,
C...  UXXXX AND UXXXXX.
C...
C...  THUS, FOR GRID POINTS 3, 4,..., N-2
C...
C...     TO RETAIN UXX
C...
C...        4*A + B + C +  4*D = 2                              (1)
C...
C...     TO DROP UXXX
C...
C...       -8*A - B + C +  8*D = 0                              (2)
C...
C...     TO DROP UXXXX
C...
C...       16*A + B + C + 16*D = 0                              (3)
C...
C...     TO DROP UXXXXX
C...
C...      -32*A - B + C + 32*D = 0                              (4)
C...
C...  EQUATIONS (1) TO (4) CAN BE SOLVED FOR A, B, C AND D.  IF EQUA-
C...  TION (1) IS ADDED TO EQUATION (2)
C...
C...        -4*A + 2*C + 12*D = 2                               (5)
C...
C...  IF EQUATION (1) IS SUBTRACTED FROM EQUATION (3)
C...
C...        12*A + 12*D = -2                                    (6)
C...
C...  IF EQUATION (1) IS ADDED TO EQUATION (4)
C...
C...        -28*A + 2*C + 36*D = 2                              (7)
C...
C...  EQUATIONS (5) TO (7) CAN BE SOLVED FOR A, C AND D.  IF EQUATION
C...  (5) IS SUBTRACTED FROM EQUATION (7), AND THE RESULT COMBINED
C...  WITH EQUATION (6)
C...
C...         12*A + 12*D = -2                                   (6)
C...
C...        -24*A + 24*D = 0                                    (8)
C...
C...  EQUATIONS (6) AND (8) CAN BE SOLVED FOR A AND D.  FROM (8), A
C...  = D.  FROM EQUATION (6), A = -1/12 AND D = -1/12.  THEN, FROM
C...  EQUATION (5), C = 4/3, AND FROM EQUATION (1), B = 4/3.
C...
C...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
C...
C...     (-1/12)*U(I-2) +   (4/3)*U(I-1) +
C...
C...       (4/3)*U(I+1) + (-1/12)*U(I+2)
C...
C...     (-1/12 + 4/3 - 1/12 + 4/3)*U(I) + UXX(I)*(DX**2) + O(DX**6)
C...
C...  OR
C...
C...     UXX(I) = (1/(12*DX**2))*(-1*U(I-2) + 16*U(I-1)
C...
C...                  - 30*U(I) + 16*U(I+1) -  1*U(I+2)         (9)
C...
C...                  + O(DX**4)
C...
C...  NOTE THAT THE UX TERM DROPS OUT, I.E., THE BASIC EQUATION IS
C...
C...        -2*A - B + C + 2*D =
C...
C...        -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0
C...
C...  EQUATION (9) WAS OBTAINED BY DROPPING ALL TERMS IN THE UNDERLYING
C...  TAYLOR SERIES UP TO AND INCLUDING THE FIFTH DERIVATIVE, UXXXXX.
C...  THUS, EQUATION (9) IS EXACT FOR POLYNOMIALS UP TO AND INCLUDING
C...  FIFTH ORDER.  THIS CAN BE CHECKED BY SUBSTITUTING THE FUNCTIONS
C...  1, X, X**2, X**3, X**4 AND X**5 IN EQUATION (9) AND COMPUTING THE
C...  CORRESPONDING DERIVATIVES FOR COMPARISON WITH THE KNOWN SECOND
C...  DERIVATIVES.  THIS IS DONE FOR 1 MERELY BY SUMMING THE WEIGHTING
C...  COEFFICIENTS IN EQUATION (9), WHICH SHOULD SUM TO ZERO, I.E.,
C...  -1 + 16 - 30 + 16 -1 = 0.
C...
C...  FOR THE REMAINING FUNCTIONS, THE ALGEBRA IS RATHER INVOLVED, BUT
C...  THESE FUNCTIONS CAN BE CHECKED NUMERICALLY, I.E., NUMERICAL VALUES
C...  OF X**2, X**3, X**4 AND X**5 CAN BE SUBSTITUTED IN EQUATION (9)
C...  AND THE COMPUTED DERIVATIVES CAN BE COMPARED WITH THE KNOW NUMERI-
C...  CAL SECOND DERIVATIVES.  THIS IS NOT A PROOF OF CORRECTNESS OF
C...  EQUATION (9), BUT WOULD LIKELY DETECT ANY ERRORS IN EQUATION (9).
C...
C...  CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (1) TO (4)
C...  CALL CHECK(1)
C...
C...  ******************************************************************
C...
C...  (2)  UXX AT THE INTERIOR POINTS I = 2 AND N-1
C...
C...  FOR GRID POINT 2, WE CONSIDER A FORMULA OF THE FORM
C...
C...     A*U(I-1) + F*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) + E*U(I+4)
C...
C...  TAYLOR SERIES EXPANSIONS OF U(I-1), U(I+1), U(I+2), U(I+3) AND
C...  U(I+4) WHEN SUBSTITUTED INTO THIS FORMULA GIVE LINEAR ALGEBRAIC
C...  EQUATIONS RELATING A, B, C, D AND E.
C...
C...     TO DROP UX
C...
C...        -A + B + 2*C + 3*D + 4*E = 0                       (10)
C...
C...     TO RETAIN UXX
C...
C...        A + B + 4*C + 9*D + 16*E = 2                       (11)
C...
C...     TO DROP UXXX
C...
C...       -A + B + 8*C + 27*D + 64*E = 0                      (12)
C...
C...     TO DROP UXXXX
C...
C...        A + B + 16*C + 81*D + 256*E = 0                    (13)
C...
C...     TO DROP UXXXXX
C...
C...       -A + B + 32*C + 243*D + 1024*E = 0                  (14)
C...
C...  EQUATIONS (11), (12), (13) AND (14) CAN BE SOLVED FOR A, B, C,
C...  D AND E.  IF EQUATION (10) IS ADDED TO EQUATION (11)
C...
C...     2*B + 6*C + 12*D +20*E = 2                            (15)
C...
C...  IF EQUATION (10) IS SUBTRACTED FROM EQUATION (12)
C...
C...     6*C + 24*D + 60*E = 0                                 (16)
C...
C...  IF EQUATION (10) IS ADDED TO EQUATION (13)
C...
C...     2*B + 18*C + 84*D + 260*E = 0                         (17)
C...
C...  IF EQUATION (10) IS SUBTRACTED FROM EQUATION (14)
C...
C...     30*C + 240*D + 1020*E = 0                             (18)
C...
C...  EQUATIONS (15), (16), (17) AND (18) CAN BE SOLVED FOR B, C, D
C...  AND E.
C...
C...     6*C + 24*D + 60*E = 0                                 (16)
C...
C...  IF EQUATION (15) IS SUBTRACTED FROM EQUATION (17)
C...
C...     12*C + 72*D + 240*E = -2                              (19)
C...
C...     30*C + 240*D + 1020*E = 0                             (18)
C...
C...  EQUATIONS (16), (18) AND (19) CAN BE SOLVED FOR C, D AND E.  IF
C...  TWO TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (19),
C...
C...     24*D + 120*E = -2                                     (20)
C...
C...  IF FIVE TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (18),
C...
C...     120*D + 720*E = 0                                     (21)
C...
C...  EQUATIONS (20) AND (21) CAN BE SOLVED FOR D AND E.  FROM (21),
C...  E = (-1/6)*D.  SUBSTITUTION IN EQUATION (20) GIVES D = -1/2.
C...  THUS, E = 1/12.  FROM EQUATION (16), C = 7/6.  FROM EQUATION
C...  (15), B = -1/3.  FROM EQUATION (10), A = 5/6.
C...
C...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
C...
C...  (5/6)*U(I-1) + (-1/3)*U(I+1) + (7/6)*U(I+2) + (-1/2)*U(I+3)
C...
C...  + (1/12)*U(I+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*U(I)
C...
C...  + UXX*(DX**2) + O(DX**6)
C...
C...  OR
C...
C...  UXX(I) = (1/12*DX**2)*(10*U(I-1) - 15*U(I) - 4*U(I+1)
C...                                                           (22)
C...         + 14*U(I+2) - 6*U(I+3) + 1*U(I+4)) + O(DX**4)
C...
C...  EQUATION (22) WILL BE APPLIED AT I = 2 AND N-1.  THUS
C...
C...  UXX(2) = (1/12*DX**2)*(10*U(1) - 15*U(2) - 4*U(3)
C...                                                           (23)
C...         + 14*U(4) - 6*U(5) + 1*U(6)) + O(DX**4)
C...
C...  UXX(N-1) = (1/12*DX**2)*(10*U(N) - 15*U(N-1) - 4*U(N-2)
C...                                                           (24)
C...           + 14*U(N-3) - 6*U(N-4) + 1*U(N-5)) + O(DX**4)
C...
C...  CHECK THE VALUES OF A, B, C, D AND E IN EQUATIONS (10) TO (14)
C...  CALL CHECK(2)
C...
C...  ******************************************************************
C...
C...  (3)  UXX AT THE BOUNDARY POINTS 1 AND N
C...
C...  FINALLY, FOR GRID POINT 1, AN APPROXIMATION WITH A NEUMANN BOUND-
C...  ARY CONDITON OF THE FORM
C...
C...     A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*UX(I) + F*U(I)
C...
C...  WILL BE USED.  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...     TO DROP UX
C...
C...        A + 2*B + 3*C + 4*D + E = 0                        (25)
C...
C...     TO RETAIN UXX
C...
C...        A + 4*B + 9*C + 16*D = 2                           (26)
C...
C...     TO DROP UXXX
C...
C...        A + 8*B + 27*C + 64*D = 0                          (27)
C...
C...     TO DROP UXXXX
C...
C...        A + 16*B + 81*C + 256*D = 0                        (28)
C...
C...     TO DROP UXXXXX
C...
C...        A + 32*B + 243*C + 1024*D = 0                      (29)
C...
C...  EQUATIONS (25) TO (29) CAN BE SOLVED FOR A, B, C, D AND E.  IF
C...
C...  EQUATION (26) IS SUBTRACTED FROM EQUATIONS (27), (28) AND (29),
C...
C...     4*B + 18*C + 48*D = -2                                (30)
C...
C...     12*B + 72*C + 240*D = -2                              (31)
C...
C...     28*B + 234*C + 1008*D = -2                            (32)
C...
C...  EQUATIONS (30), (31) AND (32) CAN BE SOLVED FOR B, C AND D
C...
C...     18*C + 96*D = 4                                       (33)
C...
C...     108*C + 672*D = 12                                    (34)
C...
C...  EQUATIONS (3) AND (34) CAN BE SOLVED FOR C AND D, C = 8/9, D =
C...  -1/8.
C...
C...  FROM EQUATION (30), B = -3.  FROM EQUATION (26), A = 8.  FROM
C...  EQUATION (25), E = -25/6.
C...
C...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
C...
C...  8*U(I+1) - 3*U(I+2) + (8/9)*U(I+3) - (1/8)*U(I+4)
C...
C...  - (25/6)*UX(I)*DX
C...
C...  = (8 - 3 + (8/9) - (1/8))*U(I) + UXX*(DX**2) + O(DX**6)
C...
C...  OR
C...
C...  UXX(I) = (1/12*DX**2)*((-415/6)*U(I) + 96*U(I+1) - 36*U(I+2)
C...                                                                (35)
C...  + (32/3)*U(I+3) - (3/2)*U(I+4) - 50*UX(I)*DX) + O(DX**4)
C...
C...  EQUATION (35) WILL BE APPLIED AT I = 1 AND I = N
C...
C...  UXX(1) = (1/12*DX**2)*((-415/6)*U(1) + 96*U(2) - 36*U(3)
C...                                                                (36)
C...  + (32/3)*U(4) - (3/2)*U(5) - 50*UX(1)*DX) + O(DX**4)
C...
C...  UXX(N) = (1/12*DX**2)*((-415/6)*U(N) + 96*U(N-1) - 36*U(N-2)
C...                                                                (37)
C...  + (32/3)*U(N-3) - (3/2)*U(N-4) + 50*UX(N)*DX) + O(DX**4)
C...
C...  CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (25) TO (29)
C...  CALL CHECK(3)
C...
C...  ALTERNATIVELY, FOR GRID POINT 1, AN APPROXIMATION WITH A DIRICHLET
C...  BOUNDARY CONDITION OF THE FORM
C...
C...  A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*U(I+5) + F*U(I)
C...
C...  CAN BE USED.  THE CORRESPONDING ALGEBRAIC EQUATIONS ARE
C...
C...     TO DROP UX
C...
C...        A + 2*B + 3*C + 4*D + 5*E = 0                      (38)
C...
C...     TO RETAIN UXX
C...
C...        A + 4*B + 9*C + 16*D + 25*E = 2                    (39)
C...
C...     TO DROP UXXX
C...
C...        A + 8*B + 27*C + 64*D + 125*E = 0                  (40)
C...
C...     TO DROP UXXXX
C...
C...        A + 16*B + 81*C + 256*D + 625*E = 0                (41)
C...
C...     TO DROP UXXXXX
C...
C...        A + 32*B + 243*C + 1024*D + 3125*E = 0             (42)
C...
C...  EQUATIONS (38), (39), (40), (41) AMD (42) CAN BE SOLVED FOR A,
C...  B, C, D AND E.
C...
C...        2*B + 6*C + 12*D + 20*E = 2                        (43)
C...
C...        6*B + 24*C + 60*D + 120*E = 0                      (44)
C...
C...        14*B + 78*C + 252*D + 620*E = 0                    (45)
C...
C...        30*B + 240*C + 1020*D + 3120*E = 0                 (46)
C...
C...  EQUATIONS (43), (44), (45) AND (46) CAN BE SOLVED FOR B, C, D
C...  AND E
C...
C...        6*C + 24*D + 60*E = -6                             (47)
C...
C...        36*C + 168*D + 480*E = -14                         (48)
C...
C...        150*C + 840*D + 2820*E = -30                       (49)
C...
C...  EQUATIONS (47), (48) AND (49) CAN BE SOLVED FOR C, D AND E
C...
C...        24*D + 120*E = 22                                  (50)
C...
C...        240*D + 1320*E = 120                               (51)
C...
C...  FROM EQUATIONS (50) AND (51), D = 61/12, E = -5/6.  FROM EQUATION
C...  (47), C = -13.  FROM EQUATION (43), B = 107/6.  FROM EQUATION
C...  (38), A = -77/6.
C...
C...  THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS
C...
C...  (-77/6)*U(I+1) + (107/6)*U(I+2) - 13*U(I+3) + (61/12)*U(I+4)
C...
C...  - (5/6)*U(I+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*U(I) +
C...
C...  UXX(I)*(DX**2) + O(DX**6)
C...
C...  OR
C...
C...  UXX(I) = (1/12*DX**2)*(45*U(I) - 154*U(I+1) + 214*U(I+2)
C...                                                                (52)
C...         - 156*U(I+3) + 61*U(I+4) - 10*U(I+5)) + O(DX**4)
C...
C...  EQUATION (52) WILL BE APPLIED AT I = 1 AND I = N
C...
C...  UXX(1) = (1/12*DX**2)*(45*U(1) - 154*U(2) + 214*U(3)
C...                                                                (53)
C...         - 156*U(4) + 61*U(5) - 10*U(6)) + O(DX**4)
C...
C...  UXX(N) = (1/12*DX**2)*(45*U(N) - 154*U(N-1) + 214*U(N-2)
C...                                                                (54)
C...         -156*U(N-3) + 61*U(N-4) - 10*U(N-5)) + O(DX**4)
C...
C...  CHECK THE VALUES OF A, B, C, D, AND E IN EQUATIONS (38) TO (42)
C...  CALL CHECK(4)
C...
C...  ******************************************************************
C...
C...  GRID SPACING
      DX=(XU-XL)/DFLOAT(N-1)
C...
C...  1/(12*DX**2) FOR SUBSEQUENT USE
      R12DXS=1./(12.0D0*DX**2)
C...
C...  UXX AT THE LEFT BOUNDARY
C...
C...     WITHOUT UX (EQUATION (53))
         IF(NL.EQ.1)THEN
         UXX(1)=R12DXS*
     1                 (   45.0D0*U(1)
     2                   -154.0D0*U(2)
     3                   +214.0D0*U(3)
     4                   -156.0D0*U(4)
     5                    +61.0D0*U(5)
     6                    -10.0D0*U(6))
C...
C...     WITH UX (EQUATION (36))
         ELSE IF(NL.EQ.2)THEN
         UXX(1)=R12DXS*
     1                 (-415.0D0/6.0D0*U(1)
     2                         +96.0D0*U(2)
     3                         -36.0D0*U(3)
     4                   +32.0D0/3.0D0*U(4)
     5                    -3.0D0/2.0D0*U(5)
     6                         -50.0D0*UX(1)*DX)
         END IF
C...
C...  UXX AT THE RIGHT BOUNDARY
C...
C...     WITHOUT UX (EQUATION (54))
         IF(NU.EQ.1)THEN
         UXX(N)=R12DXS*
     1                 (   45.0D0*U(N  )
     2                   -154.0D0*U(N-1)
     3                   +214.0D0*U(N-2)
     4                   -156.0D0*U(N-3)
     5                    +61.0D0*U(N-4)
     6                    -10.0D0*U(N-5))
C...
C...     WITH UX (EQUATION (37))
         ELSE IF(NU.EQ.2)THEN
         UXX(N)=R12DXS*
     1                 (-415.0D0/6.0D0*U(N  )
     2                         +96.0D0*U(N-1)
     3                         -36.0D0*U(N-2)
     4                   +32.0D0/3.0D0*U(N-3)
     5                    -3.0D0/2.0D0*U(N-4)
     6                         +50.0D0*UX(N  )*DX)
         END IF
C...
C...  UXX AT THE INTERIOR GRID POINTS
C...
C...     I = 2 (EQUATION (23))
         UXX(2)=R12DXS*
     1                 (   10.0D0*U(1)
     2                    -15.0D0*U(2)
     3                     -4.0D0*U(3)
     4                    +14.0D0*U(4)
     5                     -6.0D0*U(5)
     6                     +1.0D0*U(6))
C...
C...     I = N-1 (EQUATION (24))
         UXX(N-1)=R12DXS*
     1                   ( 10.0D0*U(N  )
     2                    -15.0D0*U(N-1)
     3                     -4.0D0*U(N-2)
     4                    +14.0D0*U(N-3)
     5                     -6.0D0*U(N-4)
     6                     +1.0D0*U(N-5))
C...
C...     I = 3, 4,..., N-2 (EQUATION (9))
         DO 1 I=3,N-2
         UXX(I)=R12DXS*
     1                 (   -1.0D0*U(I-2)
     2                    +16.0D0*U(I-1)
     3                    -30.0D0*U(I  )
     4                    +16.0D0*U(I+1)
     5                     -1.0D0*U(I+2))
1        CONTINUE
      RETURN
      END
