*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 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... C... SUBROUTINE DSS034 IS USED IN ESSENTIALLY THE SAME WAY AS C... SUBROUTINE DSS024. IT IS MORE GENERAL THAN DSS024 SINCE IT C... USES CENTERED OR BIASED UPWIND APPROXIMATIONS AS SELECTED BY C... THE USER. ALSO, THE CODING IN DSS034 IS MUCH SIMPLER THAN IN C... DSS024 BECAUSE CALLS TO THE ONE-DIMENSIONAL ROUTINES DSS004 C... (FIVE-POINT CENTERED) AND DSS020 (FIVE-POINT BIASED UPWIND) C... ARE USED IN PLACE OF MUCH OF THE COMPLEX CODING OF DSS024. C... BECAUSE OF THIS RELATIVE SIMPLICITY, DSS034 SHOULD BE AT LEAST C... AS EFFICIENT AS DSS024. THEREFORE, DSS034 IS RECOMMENDED FOR C... USE IN PLACE OF DSS024, 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 DERIVAT