C ALGORITHM 651, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 13, NO. 3, P. 235. SUBROUTINE HFFT3 (COEFU, PRHS, BRHS, AX, BX, AY, BY, AZ, BZ, * NX, NY, NZ, BCTY, IORDER, U, LDXU, LDYU, WORK, * NWORK, INFO) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C HFFT2 = 2 DIMENSIONS, PROBLEM DEFINED BY FUNCTIONS C HFFT2A = 2 DIMENSIONS, PROBLEM DEFINED BY ARRAYS C HFFT3 = 3 DIMENSIONS, PROBLEM DEFINED BY FUNCTIONS C HFFT3A = 3 DIMENSIONS, PROBLEM DEFINED BY ARRAYS C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C * C C P U R P O S E C ------------- C C HFFT3 SOLVES THE HELMHOLTZ EQUATION IN CARTESIAN COORDINATES C ON A THREE-DIMENSIONAL RECTANGULAR DOMAIN WITH ANY COMBINATION C OF DIRICHLET, NEUMANN, OR AND PERIODIC BOUNDARY CONDITIONS. C C C D E S C R I P T I O N C --------------------- C C HFFT3 SOLVES THE EQUATION C C 2 2 2 C D U D U D U C --- + --- + --- + COEFU*U = G C 2 2 2 C DX DY DZ C C C IN THE BOX (AX,BX)X(AY,BY)X(AZ,BZ) C C C ----------------------------. C / /: C / / : C / TOP / : C / (4) / : C / / : C / Y=BY / : C / / : C ---------------------------- : C : : : C : : : C : : RIGHT : C : : (1) : C : FRONT : : C LEFT : (5) : X=BX : C (3) : : : C : Z=BZ : / C X=AX : : / C : : / C : : / C : : / C : : / C : : / C ----------------------------/ C BOTTOM (2) Y=AY C C C WITH SOME COMBINATION OF DIRICHLET (SOLUTION PRESCRIBED), NEUMANN C (FIRST DERIVATIVE PRESCRIBED), OR PERIODIC BOUNDARY CONDITIONS. C C THE OUTPUT OF THIS PROGRAM IS A THREE-DIMENSIONAL ARRAY GIVING C ESTIMATES OF THE SOLUTION AT A SET OF GRID POINTS (X(I),Y(J),Z(K)), C FOR I=1,..,NX, J=1,..,NY, AND K=1,..,NZ, WHERE C C X(I) = AX + (I-1)*HX C Y(J) = AY + (J-1)*HY C Z(K) = AZ + (K-1)*HZ C HX = (BX-AX)/(NX-1) C HY = (BY-AY)/(NY-1) C HZ = (BZ-AZ)/(NZ-1) C C THE USER MUST CHOOSE NX, NY, AND NZ SO THAT THE GRID SPACING IS C THE SAME IN EACH OF X, Y, AND Z, I.E., HX = HY = HZ. C C WHEN COEFU=0 AND ONLY NEUMANN OR PERIODIC BOUNDARY CONDITIONS ARE C PRESCRIBED, THEN ANY CONSTANT MAY BE ADDED TO THE SOLUTION TO C OBTAIN ANOTHER SOLUTION TO THE PROBLEM. IN THIS CASE THE SOLUTION C OF MINIMUM INFINITY NORM IS RETURNED. C C THE SOLUTION IS COMPUTED USING A FOURTH ORDER ACCURATE FINITE C DIFFERENCE APPROXIMATION OF THE CONTINUOUS EQUATION. THE RESULTING C SYSTEM OF LINEAR ALGEBRAIC EQUATIONS IS SOLVED USING FAST FOURIER C TRANSFORM TECHNIQUES. THE ALGORITHM RELIES UPON THE FACT THAT C NX-1 AND NZ-1 ARE HIGHLY COMPOSITE (THE PRODUCT OF SMALL PRIMES). C C C P A R A M E T E R S C ------------------- C C COEFU REAL SCALAR (INPUT) C THE COEFFICIENT OF U IN THE PARTIAL DIFFERENTIAL EQUATION. C C PRHS REAL FUNCTION OF X,Y,Z (INPUT) C RETURNS THE RIGHT-HAND SIDE OF THE DIFFERENTIAL EQUATION C FOR ANY (X,Y,Z) IN THE DOMAIN OR ON THE BOUNDARY (SEE C DESCRIPTION BELOW). THE NAME OF THIS FUNCTION MUST BE C DECLARED EXTERNAL IN THE CALLING PROGRAM. C C BRHS REAL FUNCTION OF K,X,Y,Z (INPUT) C RETURNS THE RIGHT-HAND SIDE OF THE BOUNDARY CONDITION AT C THE POINT (X,Y,Z) ON THE K-TH SIDE OF THE DOMAIN (SEE C DESCRIPTION BELOW). THE NAME OF THIS FUNCTION MUST BE C DECLARED EXTERNAL IN THE CALLING PROGRAM. C C AX REAL SCALAR (INPUT) .LT. BX C THE VALUE OF X ALONG THE LEFT SIDE OF THE DOMAIN. C C BX REAL SCALAR (INPUT) C THE VALUE OF X ALONG THE RIGHT SIDE OF THE DOMAIN. C C AY REAL SCALAR (INPUT) .LT. BY C THE VALUE OF Y ALONG THE BOTTOM SIDE OF THE DOMAIN. C C BY REAL SCALAR (INPUT) C THE VALUE OF Y ALONG THE TOP SIDE OF THE DOMAIN. C C AZ REAL SCALAR (INPUT) .LT. BZ C THE VALUE OF Z ALONG THE FRONT SIDE OF THE DOMAIN. C C BZ REAL SCALAR (INPUT) C THE VALUE OF Z ALONG THE BACK SIDE OF THE DOMAIN. C C NX INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN X. C (THE NUMBER OF SUB-INTERVALS IN X IS THEN NX-1.) C NX SHOULD BE CHOSEN SO THAT NX-1 IS HIGHLY COMPOSITE C (THE PRODUCT OF SMALL PRIMES) TO INCREASE THE SPEED C OF THE FOURIER TRANSFORM ALGORITHM. C C NY INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN Y. C (THE NUMBER OF SUB-INTERVALS IN Y IS THEN NY-1.) C C NZ INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN Z. C (THE NUMBER OF SUB-INTERVALS IN Z IS THEN NZ-1.) C NZ SHOULD BE CHOSEN SO THAT NX-1 IS HIGHLY COMPOSITE C (THE PRODUCT OF SMALL PRIMES) TO INCREASE THE SPEED C OF THE FOURIER TRANSFORM ALGORITHM. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C INDICATES TYPE OF BOUNDARY CONDITION ON EACH SIDE OF THE C DOMAIN AS FOLLOWS. C C BCTY(1) = TYPE ON RIGHT SIDE (X=BX) C BCTY(2) = TYPE ON BOTTOM SIDE (Y=AY) C BCTY(3) = TYPE ON LEFT SIDE (X=AX) C BCTY(4) = TYPE ON TOP SIDE (Y=BY) C BCTY(5) = TYPE ON FRONT SIDE (Z=BZ) C BCTY(6) = TYPE ON BACK SIDE (Z=AZ) C C POSSIBLE VALUES ARE C C 1 == U PRESCRIBED C 2 == DU/DX PRESCRIBED (FOR BCTY(1) AND BCTY(3)) OR C DU/DY PRESCRIBED (FOR BCTY(2) AND BCTY(4)) OR C DU/DZ PRESCRIBED (FOR BCTY(5) AND BCTY(6)) C 3 == PERIODIC C C IORDER INTEGER SCALAR (INPUT) C THE ORDER OF ACCURACY OF THE FINITE DIFFERENCE C APPROXIMATION TO THE PROBLEM. POSSIBLE VALUES ARE C C 2 == SECOND ORDER ACCURATE COMPACT 9-POINT DIFFERENCES C 4 == FOURTH ORDER ACCURATE COMPACT 9-POINT DIFFERENCES C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ+2 (INPUT/OUTPUT) C ON INPUT, U(I,J,K) CONTAINS THE VALUE OF THE RIGHT HAND C SIDE OF THE PARTIAL DIFFERENTAIL EQUATION AT THE (I,J,K)TH C GRID POINT. THAT IS, C C U(I,J,K) = G(X(I),Y(J),Z(K)) C C FOR I=1,..,NX, J=1,..,NY, AND K=1,..,NZ. C ON OUTPUT, U(I,J,K) IS THE VALUE OF THE COMPUTED SOLUTION C AT THE POINT (X(I),Y(J),Z(K)). THE PLANES U(NX+1,*,*), C U(NX+2,*,*), U(*,NY+1,*), U(*,NY+2,*), U(*,*,NZ+1), AND C U(*,*,NZ+2) ARE USED AS WORKING STORAGE. C C LDXU INTEGER SCALAR (INPUT) .GE. NX+2 C THE LEADING DIMENSION OF THE ARRAY U EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C LDYU INTEGER SCALAR (INPUT) .GE. NY+2 C THE SECOND DIMENSION OF THE ARRAY U EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C WORK REAL ARRAY OF SIZE NWORK. C WORKING STORAGE REQUIRED BY HFFT3. C C NWORK INTEGER SCALAR (INPUT) C .GE. (NX+1)*(NY+1)*(NZ+1)*(IORDER-2)/2 C + 2*(NX*NY+NX*NZ+NY*NZ) C + MAX( (NX+1)*(NY+1)*(IORDER-2), C (NX+3)*(NZ+5) + 5*NY + (NX+NZ)/2 + 15 ) C THE LENGTH OF THE ARRAY WORK AS DECLARED IN THE CALLING C PROGRAM. THIS MAY BE REDUCED BY 4*NY IF COEFU.LE.0. C C INFO INTEGER SCALAR (OUTPUT) C INDICATES STATUS OF COMPUTED SOLUTION. C POSSIBLE VALUES ARE C C 2 == WARNING. NO SOLUTION EXISTS UNLESS A CONSISTENCY C CONDITION IS SATISFIED (SEE REF. 4). THE DISCRETE C PROBLEM IS ADJUSTED (BY ADDING A CONSTANT TO THE C RIGHT SIDE) SO THAT THIS CONDITION IS SATISFIED. C THIS CONSTANT IS RETURNED IN WORK(1). IF IT IS NOT C SMALL THEN THE PROBLEM MAY NOT BE WELL-POSED. IN C ADDITION, THE SOLUTION IS UNIQUE ONLY UP TO AN C ADDITIVE CONSTANT. C 1 == WARNING. COEFU.GT.0 A SOLUTION MAY NOT EXIST IF C COEFU IS AN EIGENVALUE OF THE LAPLACIAN. IF COEFU C IS NEAR ONE OF THESE VALUES THEN THE COMPUTED C SOLUTION MAY BE UNRELIABLE. C 0 == SUCCESS. SUBPROGRAM RAN TO COMPLETION. C -1 == ERROR. NX.LT.4 C -2 == ERROR. NY.LT.4 C -3 == ERROR. NZ.LT.4 C -4 == ERROR. LDXU.LT.NX+2 C -5 == ERROR. LDYU.LT.NY+2 C -6 == ERROR. IORDER NOT 2 OR 4. C -7 == ERROR. ELEMENT OF BCTY NOT 1, 2 OR 3. C -8 == ERROR. PERIODIC BOUNDARY CONDITIONS SPECIFIED ON C ONE SIDE OF DOMAIN BUT NOT ON THE OPPOSITE. C -9 == ERROR. NWORK TOO SMALL. C -10 == ERROR. BX.LT.AX C -11 == ERROR. BY.LT.AY C -12 == ERROR. BZ.LT.AZ C -13 == ERROR. GRID SIZE IN Y NOT SAME AS IN X. C -14 == ERROR. GRID SIZE IN Z NOT SAME AS IN X. C C C U S E R - D E F I N E D F U N C T I O N S C ------------------------------------------- C C THE RIGHT HAND SIDE (FORCING FUNCTION) OF THE PARTIAL C DIFFERENTIAL EQUATION IS SPECIFIED BY THE FUNCTION PRHS. C C REAL FUNCTION PRHS (X,Y,Z) C REAL X,Y,Z C PRHS = RIGHT SIDE OF THE PDE AT (X,Y,Z) C RETURN C END C C THE VALUE OF THE SOLUTION OR ITS FIRST DERIVATIVE ALONG THE C EDGES OF THE DOMAIN WHERE BCTY=1 OR 2 IS SUPPLIED VIA THE C FUNCTION BRHS. C C REAL FUNCTION BRHS (K,X,Y,Z) C INTEGER K C REAL X,Y,Z C BRHS = RIGHT HAND SIDE OF KTH BOUNDARY CONDITION C (SEE BCTY) AT THE POINT (X,Y,Z) C RETURN C END C C C E X T E R N A L R E F E R E N C E S C ------------------------------------- C C HFFT3A,FDIS3,HDIS3,FD3N, C STORD3,FD3D,FD2D,FD2DA,HD3N, C REFL3,MDALG3,MDALG2,EVDISC, C TRISOL,TRSALL,TRSOLG,TRSOLP, C SGPSL,FFTI,FFTB,FFTF -- THIS PACKAGE C C RFFTI,RFFTF,RFFTB, C SINTI,SINT,COSTI,COST, C SINQI,SINQF,SINQB, C COSQI,COSQF,COSQB -- FFTPACK (SEE REF. 2) C C R1MACH -- MACHINE CONSTANTS (SEE REF. 3) C C PRHS,BRHS -- USER-SUPPLIED C C C P O R T A B I L I T Y C --------------------- C C THIS PACKAGE IS WRITTEN IN ANSI STANDARD FORTRAN (1977). C ALL MACHINE-DEPENDENT QUANTITIES ARE OBTAINED FROM THE C FUNCTION R1MACH (SEE REF. 3). C C C R E F E R E N C E S C ------------------- C C 1) R. BOISVERT, A FOURTH ORDER ACCURATE FAST DIRECT METHOD C FOR THE HELMHOLTZ EQUATION, IN ELLIPTIC PROBLEM SOLVERS C II (G. BIRKHOFF AND A. SCHOENSTADT, EDS.), ACADEMIC PRESS, C ORLANDO, FLA., 1984, 35-44. C C 2) THE FFT PACKAGE USED IS A SLIGHTLY MODIFIED VERSION OF THE C PACKAGE FFTPACK WRITTEN BY PAUL SWARZTRAUBER. FFTPACK IS C ALSO USED BY THE SUBPROGRAM HW3CRT OF FISHPAK (VERSION 3). C FFTPACK IS DESCRIBED IN C C P.N. SWARZTRAUBER, VECTORIZING THE FFTS, IN PARALLEL C COMPUTATIONS (G. RODRIGUE, ED.), ACADEMIC PRESS, 1982, C PP. 51-83. C C FOR FURTHER INFORMATION WRITE INFORMATION SERVICES OFFICE, C COMPUTING FACILITY, NATIONAL CENTER FOR ATMOSPHERIC RESEARCH, C BOX 3000, BOULDER, CO 80303, USA. C C 3) P. FOX, A. HALL, AND N. SCHRYER, ALGORITHM 528: FRAMEWORK C FOR A PORTABLE LIBRARY, ACM TRANS. MATH. SOFT. 4 (1978), C PP. 177-188. C C 4) S. G. MIKHLIN (ED.), LINEAR EQUATIONS OF MATHEMATICAL PHYSICS, C HOLT, RINEHART AND WINSTON, NEW YORK, 1967. C C C A U T H O R / V E R S I O N C ----------------------------- C C RONALD F. BOISVERT C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C GAITHERSBURG, MD 20899 C USA C C ORIGINAL DECEMBER 1985 C REVISED APRIL 1987 C C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER NX, NY, NZ, BCTY(6), IORDER, LDXU, LDYU, NWORK, INFO REAL * COEFU, AX, BX, AY, BY, AZ, BZ, U(LDXU,LDYU,*), WORK(NWORK) C C ... LOCAL VARIABLES C INTEGER LOCGH, LOCBD1, LOCBD2, LOCBD3, LOCBD4, LOCBD5, LOCBD6, * NEEDED, NWORKA, NXM1, NYM1, NZM1, NXP1, NYP1, NZP1 REAL * EPSM, R1MACH, PRHS, BRHS, H, HX, HY, HZ, AXH, AYH, AZH, * X, Y, Z C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C EPSM = R1MACH(4) C NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 NXP1 = NX + 1 NYP1 = NY + 1 NZP1 = NZ + 1 C H = (BX-AX)/REAL(NXM1) C LOC GH = 1 LOC BD1 = LOC GH IF (IORDER .EQ. 4) LOC BD1 = LOC GH + NXP1*NYP1*NZP1 LOC BD2 = LOC BD1 + NY*NZ LOC BD3 = LOC BD2 + NX*NZ LOC BD4 = LOC BD3 + NY*NZ LOC BD5 = LOC BD4 + NX*NZ LOC BD6 = LOC BD5 + NX*NY LOC WRK = LOC BD6 + NX*NY NEEDED = LOC WRK + MAX( NXP1*NXP1*(IORDER-2), * (NX+3)*(NY+5) + 5*NY + (NX+NZ)/2 + 15 ) - 1 NWORKA = NWORK - LOCWRK + 1 C C C --------------------------- C CHECK VALIDITY OF ARGUMENTS C --------------------------- C INFO = 0 IF (NX .LT. PRDC) GO TO 901 IF (NY .LT. PRDC) GO TO 902 IF (NZ .LT. PRDC) GO TO 903 IF (LDXU .LT. NX+2) GO TO 904 IF (LDYU .LT. NY+2) GO TO 905 IF ((IORDER .NE. 2) .AND. (IORDER .NE. 4)) GO TO 906 DO 10 K=1,6 IF ((BCTY(K) .LT. DRCH) .OR. (BCTY(K) .GT. PRDC)) GO TO 907 10 CONTINUE IF (((BCTY(RIGHT ) .EQ. PRDC) .AND. (BCTY(LEFT ) .NE. PRDC)) .OR. * ((BCTY(LEFT ) .EQ. PRDC) .AND. (BCTY(RIGHT ) .NE. PRDC)) .OR. * ((BCTY(BOTTOM) .EQ. PRDC) .AND. (BCTY(TOP ) .NE. PRDC)) .OR. * ((BCTY(TOP ) .EQ. PRDC) .AND. (BCTY(BOTTOM) .NE. PRDC)) .OR. * ((BCTY(FRONT ) .EQ. PRDC) .AND. (BCTY(BACK ) .NE. PRDC)) .OR. * ((BCTY(BACK ) .EQ. PRDC) .AND. (BCTY(FRONT ) .NE. PRDC)) ) * GO TO 908 IF (NWORK .LT. NEEDED) GO TO 909 IF (BX .LE. AX) GO TO 910 IF (BY .LE. AY) GO TO 911 IF (BZ .LE. AZ) GO TO 912 HX = H HY = (BY-AY)/REAL(NYM1) HZ = (BZ-AZ)/REAL(NZM1) IF (ABS(HX-HY) .GT. EPSM) GO TO 913 IF (ABS(HX-HZ) .GT. EPSM) GO TO 914 C C C ------------------ C COMPUTE RHS OF PDE C ------------------ C C ... AT GRID POINTS C DO 100 K=1,NZ Z = AZ + REAL(K-1)*H DO 100 J=1,NY Y = AY + REAL(J-1)*H DO 100 I=1,NX X = AX + REAL(I-1)*H U(I,J,K) = PRHS(X,Y,Z) 100 CONTINUE C C ... AT HALF GRID POINTS C IF (IORDER .EQ. 4) THEN AXH = AX + 0.50E0*H AYH = AY + 0.50E0*H AZH = AZ + 0.50E0*H DO 200 K=1,NZM1 Z = AZH + REAL(K-1)*H DO 200 J=1,NYM1 Y = AYH + REAL(J-1)*H LOC = LOC GH + (K*NYP1 + J)*NXP1 DO 200 I=1,NXM1 X = AXH + REAL(I-1)*H LOC = LOC + 1 WORK(LOC) = PRHS(X,Y,Z) 200 CONTINUE ENDIF C C C ------------------- C STORE BOUNDARY DATA C ------------------- C C ... BOTTOM SIDE C IF (BCTY(BOTTOM) .NE. PRDC) THEN DO 265 K=1,NZ Z = AZ + REAL(K-1)*H LOC = LOC BD2 + (K-1)*NX - 1 DO 265 I=1,NX X = AX + REAL(I-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(BOTTOM,X,AY,Z) 265 CONTINUE ENDIF C C ... TOP SIDE C IF (BCTY(TOP) .NE. PRDC) THEN DO 275 K=1,NZ Z = AZ + REAL(K-1)*H LOC = LOC BD4 + (K-1)*NX - 1 DO 275 I=1,NX X = AX + REAL(I-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(TOP,X,BY,Z) 275 CONTINUE ENDIF C C ... LEFT SIDE C IF (BCTY(LEFT) .NE. PRDC) THEN DO 285 K=1,NZ Z = AZ + REAL(K-1)*H LOC = LOC BD3 + (K-1)*NY - 1 DO 285 J=1,NY Y = AY + REAL(J-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(LEFT,AX,Y,Z) 285 CONTINUE ENDIF C C ... RIGHT SIDE C IF (BCTY(RIGHT) .NE. PRDC) THEN DO 295 K=1,NZ Z = AZ + REAL(K-1)*H LOC = LOC BD1 + (K-1)*NY - 1 DO 295 J=1,NY Y = AY + REAL(J-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(RIGHT,BX,Y,Z) 295 CONTINUE ENDIF C C ... FRONT SIDE C IF (BCTY(FRONT) .NE. PRDC) THEN DO 305 J=1,NY Y = AY + REAL(J-1)*H LOC = LOC BD5 + (J-1)*NX - 1 DO 305 I=1,NX X = AX + REAL(I-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(FRONT,X,Y,BZ) 305 CONTINUE ENDIF C C ... BACK SIDE C IF (BCTY(BACK) .NE. PRDC) THEN DO 315 J=1,NY Y = AY + REAL(J-1)*H LOC = LOC BD6 + (J-1)*NX - 1 DO 315 I=1,NX X = AX + REAL(I-1)*H LOC = LOC + 1 WORK(LOC) = BRHS(BACK,X,Y,AZ) 315 CONTINUE ENDIF C C C ---------------- C CALL FAST SOLVER C ---------------- C CALL HFFT3A(COEFU,NX,NY,NZ,H,WORK(LOCGH),NXP1,NYP1,BCTY, * WORK(LOCBD1),WORK(LOCBD2),WORK(LOCBD3),WORK(LOCBD4), * WORK(LOCBD5),WORK(LOCBD6),NX,NY,IORDER,U,LDXU,LDYU, * WORK(LOCWRK),NWORKA,INFO) C C C ----------- C NORMAL EXIT C ----------- C WORK(1) = WORK(LOCWRK) GO TO 999 C C C ----------- C ERROR EXITS C ----------- C 901 CONTINUE INFO = -1 GO TO 999 C 902 CONTINUE INFO = -2 GO TO 999 C 903 CONTINUE INFO = -3 GO TO 999 C 904 CONTINUE INFO = -4 GO TO 999 C 905 CONTINUE INFO = -5 GO TO 999 C 906 CONTINUE INFO = -6 GO TO 999 C 907 CONTINUE INFO = -7 GO TO 999 C 908 CONTINUE INFO = -8 GO TO 999 C 909 CONTINUE INFO = -9 GO TO 999 C 910 CONTINUE INFO = -10 GO TO 999 C 911 CONTINUE INFO = -11 GO TO 999 C 912 CONTINUE INFO = -12 GO TO 999 C 913 CONTINUE INFO = -13 GO TO 999 C 914 CONTINUE INFO = -14 GO TO 999 C 999 CONTINUE RETURN END SUBROUTINE HFFT3A (COEFU, NX, NY, NZ, H, GH, LDXGH, LDYGH, BCTY, * BD1, BD2, BD3, BD4, BD5, BD6, LDXBD, LDYBD, * IORDER, U, LDXU, LDYU, WORK, NWORK, INFO) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C HFFT2 = 2 DIMENSIONS, PROBLEM DEFINED BY FUNCTIONS C HFFT2A = 2 DIMENSIONS, PROBLEM DEFINED BY ARRAYS C HFFT3 = 3 DIMENSIONS, PROBLEM DEFINED BY FUNCTIONS C HFFT3A = 3 DIMENSIONS, PROBLEM DEFINED BY ARRAYS C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C P U R P O S E C ------------- C C HFFT3A SOLVES THE HELMHOLTZ EQUATION IN CARTESIAN COORDINATES C ON A THREE-DIMENSIONAL RECTANGULAR DOMAIN WITH ANY COMBINATION C OF DIRICHLET, NEUMANN, OR AND PERIODIC BOUNDARY CONDITIONS. C C C D E S C R I P T I O N C --------------------- C C HFFT3A SOLVES THE EQUATION C C 2 2 2 C D U D U D U C --- + --- + --- + COEFU*U = G C 2 2 2 C DX DY DZ C C C IN THE BOX (AX,BX)X(AY,BY)X(AZ,BZ) C C C ----------------------------. C / /: C / / : C / TOP / : C / (4) / : C / / : C / Y=BY / : C / / : C ---------------------------- : C : : : C : : : C : : RIGHT : C : : (1) : C : FRONT : : C LEFT : (5) : X=BX : C (3) : : : C : Z=BZ : / C X=AX : : / C : : / C : : / C : : / C : : / C : : / C ---------------------------- C BOTTOM (2) Y=AY C C C WITH SOME COMBINATION OF DIRICHLET (SOLUTION PRESCRIBED), NEUMANN C (FIRST DERIVATIVE PRESCRIBED), OR PERIODIC BOUNDARY CONDITIONS. C C THE OUTPUT OF THIS PROGRAM IS A THREE-DIMENSIONAL ARRAY GIVING C ESTIMATES OF THE SOLUTION AT A SET OF GRID POINTS (X(I),Y(J),Z(K)), C FOR I=1,..,NX, J=1,..,NY, AND K=1,..,NZ, WHERE C C X(I) = AX + (I-1)*H C Y(J) = AY + (J-1)*H C Z(K) = AZ + (K-1)*H C C WHEN COEFU=0 AND ONLY NEUMANN OR PERIODIC BOUNDARY CONDITIONS ARE C PRESCRIBED, THEN ANY CONSTANT MAY BE ADDED TO THE SOLUTION TO C OBTAIN ANOTHER SOLUTION TO THE PROBLEM. IN THIS CASE THE SOLUTION C OF MINIMUM INFINITY NORM IS RETURNED. C C THE SOLUTION IS COMPUTED USING A FOURTH ORDER ACCURATE FINITE C DIFFERENCE APPROXIMATION OF THE CONTINUOUS EQUATION. THE RESULTING C SYSTEM OF LINEAR ALGEBRAIC EQUATIONS IS SOLVED USING FAST FOURIER C TRANSFORM TECHNIQUES. THE ALGORITHM RELIES UPON THE FACT THAT C NX-1 AND NZ-1 ARE HIGHLY COMPOSITE (THE PRODUCT OF SMALL PRIMES). C C C P A R A M E T E R S C ------------------- C C COEFU REAL SCALAR (INPUT) C THE COEFFICIENT OF U IN THE PARTIAL DIFFERENTIAL EQUATION. C C NX INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN X. C (THE NUMBER OF SUB-INTERVALS IN X IS THEN NX-1.) C NX SHOULD BE CHOSEN SO THAT NX-1 IS HIGHLY COMPOSITE C (THE PRODUCT OF SMALL PRIMES) TO INCREASE THE SPEED C OF THE FOURIER TRANSFORM ALGORITHM. C C NY INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN Y. C (THE NUMBER OF SUB-INTERVALS IN Y IS THEN NY-1.) C C NZ INTEGER SCALAR (INPUT) .GE. 4 C THE NUMBER OF GRID LINES IN Z. C (THE NUMBER OF SUB-INTERVALS IN Z IS THEN NZ-1.) C NZ SHOULD BE CHOSEN SO THAT NX-1 IS HIGHLY COMPOSITE C (THE PRODUCT OF SMALL PRIMES) TO INCREASE THE SPEED C OF THE FOURIER TRANSFORM ALGORITHM. C C H REAL SCALAR (INPUT) C THE SPACE BETWEEN ADJACENT GRID LINES (IN X, Y, AND Z). C C GH REAL ARRAY OF SIZE LDXGH BY LDYGH BY NZ+1 (INPUT) C THE RIGHT HAND SIDE OF THE PDE AT HALF GRID POINTS. C G(I+1,J+1,K+1) IS THE VALUE OF THE RIGHT HAND SIDE AT THE C (I,J,K)TH HALF GRID POINT, I=1,..,NX-1, J=1,..,NY-1, C K=1,..,NZ-1. THAT IS, C C GH(I+1,J+1,K+1) = G(XH(I),YH(J),ZH(K)) C XH(I) = AX + (I-0.5)*H C YH(J) = AY + (J-0.5)*H C ZH(K) = AZ + (K-0.5)*H C C PLANES GH(1,*,*), GH(NX+1,*,*), GH(*,1,*), GH(*,NY+1,*), C GH(*,*,1), AND GH(*,*,NZ+2) ARE USED FOR WORKING STORAGE. C GH IS NOT USED WHEN IORDER=2. C C LDXGH INTEGER SCALAR (INPUT) .GE. NX+1 C THE LEADING DIMENSION OF THE ARRAY GH EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C LDYGH INTEGER SCALAR (INPUT) .GE. NY+1 C THE SECOND DIMENSION OF THE ARRAY GH EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C INDICATES TYPE OF BOUNDARY CONDITION ON EACH SIDE OF THE C DOMAIN AS FOLLOWS. C C BCTY(1) = TYPE ON RIGHT SIDE (X=BX) C BCTY(2) = TYPE ON BOTTOM SIDE (Y=AY) C BCTY(3) = TYPE ON LEFT SIDE (X=AX) C BCTY(4) = TYPE ON TOP SIDE (Y=BY) C BCTY(5) = TYPE ON FRONT SIDE (Z=BZ) C BCTY(6) = TYPE ON BACK SIDE (Z=AZ) C C POSSIBLE VALUES ARE C C 1 == U PRESCRIBED C 2 == DU/DX PRESCRIBED (FOR BCTY(1) AND BCTY(3)) OR C DU/DY PRESCRIBED (FOR BCTY(2) AND BCTY(4)) OR C DU/DZ PRESCRIBED (FOR BCTY(5) AND BCTY(6)) C 3 == PERIODIC C C BD1 REAL ARRAY OF SIZE LDYBD BY NZ (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE RIGHT SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(1) VALUE STORED C ------------------------------ C 1 U C 2 DU/DX C 3 NONE C C BD1(J,K) GIVES THE VALUE AT THE POINT (BX,Y(J),Z(K)). C C BD2 REAL ARRAY OF SIZE LDXBD BY NZ (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE BOTTOM SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(2) VALUE STORED C ------------------------------ C 1 U C 2 DU/DY C 3 NONE C C BD2(I,K) GIVES THE VALUE AT THE POINT (X(I),AY,Z(K)). C C BD3 REAL ARRAY OF SIZE LDYBD BY NZ (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE LEFT SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(3) VALUE STORED C ------------------------------ C 1 U C 2 DU/DX C 3 NONE C C BD3(J,K) GIVES THE VALUE AT THE POINT (AX,Y(J),Z(K)). C C BD4 REAL ARRAY OF SIZE LDXBD BY NZ (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE TOP SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(4) VALUE STORED C ------------------------------ C 1 U C 2 DU/DY C 3 NONE C C BD4(I,K) GIVES THE VALUE AT THE POINT (X(I),BY,Z(K)). C C BD5 REAL ARRAY OF SIZE LDXBD BY NY (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE FRONT SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(5) VALUE STORED C ------------------------------ C 1 U C 2 DU/DZ C 3 NONE C C BD5(I,J) GIVES THE VALUE AT THE POINT (X(I),Y(J),BZ). C C BD6 REAL ARRAY OF SIYE LDXBD BY NY (INPUT) C THE VALUES OF THE BOUNDARY CONDITION AT GRID POINTS ON C THE BACK SIDE OF THE DOMAIN. THE VALUE STORED DEPENDS C UPON THE TYPE OF BOUNDARY CONDITION AS FOLLOWS C C BCTY(6) VALUE STORED C ------------------------------ C 1 U C 2 DU/DZ C 3 NONE C C BD6(I,J) GIVES THE VALUE AT THE POINT (X(I),Y(J),AZ). C C LDXBD INTEGER SCALAR (INPUT) .GE. NX C THE LEADING DIMENSION OF THE ARRAYS BD2, BD4, BD5, AND C BD6, EXACTLY AS SPECIFIED IN THE CALLING PROGRAM. C C LDYBD INTEGER SCALAR (INPUT) .GE. NY C THE LEADING DIMENSION OF THE ARRAYS BD1 AND BD3, EXACTLY C AS SPECIFIED IN THE CALLING PROGRAM. C C IORDER INTEGER SCALAR (INPUT) C THE ORDER OF ACCURACY OF THE FINITE DIFFERENCE C APPROXIMATION TO THE PROBLEM. POSSIBLE VALUES ARE C C 2 == SECOND ORDER ACCURATE 9-POINT COMPACT DIFFERENCES C 4 == FOURTH ORDER ACCURATE 9-POINT COMPACT DIFFERENCES C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ+2 (INPUT/OUTPUT) C ON INPUT, U(I,J,K) CONTAINS THE VALUE OF THE RIGHT HAND C SIDE OF THE PARTIAL DIFFERENTAIL EQUATION AT THE (I,J,K)TH C GRID POINT. THAT IS, C C U(I,J,K) = G(X(I),Y(J),Z(K)) C C FOR I=1,..,NX, J=1,..,NY, AND K=1,..,NZ. C ON OUTPUT, U(I,J,K) IS THE VALUE OF THE COMPUTED SOLUTION C AT THE POINT (X(I),Y(J),Z(K)). THE PLANES U(NX+1,*,*), C U(NX+2,*,*), U(*,NY+1,*), U(*,NY+2,*), U(*,*,NZ+1), AND C U(*,*,NZ+2) ARE USED AS WORKING STORAGE. C C LDXU INTEGER SCALAR (INPUT) .GE. NX+2 C THE LEADING DIMENSION OF THE ARRAY U EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C LDYU INTEGER SCALAR (INPUT) .GE. NY+2 C THE SECOND DIMENSION OF THE ARRAY U EXACTLY AS SPECIFIED C IN THE CALLING PROGRAM. C C WORK REAL ARRAY OF SIZE NWORK. C WORKING STORAGE REQUIRED BY HFFT3A. C C NWORK INTEGER SCALAR (INPUT) C .GE. MAX( (NX+1)*(NY+1)*(IORDER-2), C (NX+3)*(NZ+5) + 5*NY + (NX+NZ)/2 + 15 ) C THE LENGTH OF THE ARRAY WORK AS DECLARED IN THE CALLING C PROGRAM. THIS MAY BE REDUCED BY 4*NY IF COEFU.LE.0. C C INFO INTEGER SCALAR (OUTPUT) C INDICATES STATUS OF COMPUTED SOLUTION. C POSSIBLE VALUES ARE C C 2 == WARNING. NO SOLUTION EXISTS UNLESS A CONSISTENCY C CONDITION IS SATISFIED (SEE REF. 4). THE DISCRETE C PROBLEM IS ADJUSTED (BY ADDING A CONSTANT TO THE C RIGHT SIDE) SO THAT THIS CONDITION IS SATISFIED. C THIS CONSTANT IS RETURNED IN WORK(1). IF IT IS NOT C SMALL THEN THE PROBLEM MAY NOT BE WELL-POSED. IN C ADDITION, THE SOLUTION IS UNIQUE ONLY UP TO AN C ADDITIVE CONSTANT. C 1 == WARNING. COEFU.GT.0 A SOLUTION MAY NOT EXIST IF C COEFU IS AN EIGENVALUE OF THE LAPLACIAN. IF COEFU C IS NEAR ONE OF THESE VALUES THEN THE COMPUTED C SOLUTION MAY BE UNRELIABLE. C 0 == SUCCESS. SUBPROGRAM RAN TO COMPLETION. C -1 == ERROR. NX.LT.4 C -2 == ERROR. NY.LT.4 C -3 == ERROR. NZ.LT.4 C -4 == ERROR. LDXU.LT.NX+2 C -5 == ERROR. LDYU.LT.NY+2 C -6 == ERROR. IORDER NOT 2 OR 4. C -7 == ERROR. ELEMENT OF BCTY NOT 1, 2 OR 3. C -8 == ERROR. PERIODIC BOUNDARY CONDITIONS SPECIFIED ON C ONE SIDE OF DOMAIN BUT NOT ON THE OPPOSITE. C -9 == ERROR. NWORK TOO SMALL. C -10 == ERROR. LDXGH.LT.NX+1 C -11 == ERROR. LDYGH.LT.NY+1 C C C C E X T E R N A L R E F E R E N C E S C ------------------------------------- C C FDIS3,HDIS3,FD3N,STORD3, C FD3D,FD2D,FD2DA,HD3N,REFL3, C MDALG3,MDALG2,EVDISC,TRISOL, C TRSALL,TRSOLG,TRSOLP, C SGPSL,FFTI,FFTB,FFTF -- THIS PACKAGE C C RFFTI,RFFTF,RFFTB, C SINTI,SINT,COSTI,COST, C SINQI,SINQF,SINQB, C COSQI,COSQF,COSQB -- FFTPACK (SEE REF. 2) C C R1MACH -- MACHINE CONSTANTS (SEE REF. 3) C C C P O R T A B I L I T Y C --------------------- C C THIS PACKAGE IS WRITTEN IN ANSI STANDARD FORTRAN (1977). C ALL MACHINE-DEPENDENT QUANTITIES ARE OBTAINED FROM THE C FUNCTION R1MACH (SEE REF. 3). C C C R E F E R E N C E S C ------------------- C C 1) R. BOISVERT, A FOURTH ORDER ACCURATE FAST DIRECT METHOD C FOR THE HELMHOLTZ EQUATION, IN ELLIPTIC PROBLEM SOLVERS C II (G. BIRKHOFF AND A. SCHOENSTADT, EDS.), ACADEMIC PRESS, C ORLANDO, FLA., 1984, 35-44. C C 2) THE FFT PACKAGE USED IS A SLIGHTLY MODIFIED VERSION OF THE C PACKAGE FFTPACK WRITTEN BY PAUL SWARZTRAUBER. FFTPACK IS C ALSO USED BY THE SUBPROGRAM HW3CRT OF FISHPAK (VERSION 3). C FFTPACK IS DESCRIBED IN C C P.N. SWARZTRAUBER, VECTORIZING THE FFTS, IN PARALLEL C COMPUTATIONS (G. RODRIGUE, ED.), ACADEMIC PRESS, 1982, C PP. 51-83. C C FOR FURTHER INFORMATION WRITE INFORMATION SERVICES OFFICE, C COMPUTING FACILITY, NATIONAL CENTER FOR ATMOSPHERIC RESEARCH, C BOX 3000, BOULDER, CO 80303, USA. C C 3) P. FOX, A. HALL, AND N. SCHRYER, ALGORITHM 528: FRAMEWORK C FOR A PORTABLE LIBRARY, ACM TRANS. MATH. SOFT. 4 (1978), C PP. 177-188. C C 4) S. G. MIKHLIN (ED.), LINEAR EQUATIONS OF MATHEMATICAL PHYSICS, C HOLT, RINEHART AND WINSTON, NEW YORK, 1967. C C C A U T H O R / V E R S I O N C ----------------------------- C C RONALD F. BOISVERT C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C GAITHERSBURG, MD 20899 C USA C C ORIGINAL DECEMBER 1985 C REVISED APRIL 1987 C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER NX, NY, NZ, LDXGH, LDYGH, LDXBD, LDYBD, IORDER, * LDXU, LDYU, BCTY(6) REAL * H, COEFU, GH(LDXGH,LDYGH,*), U(LDXU,LDYU,*), * BD1(LDYBD,NZ), BD2(LDXBD,NZ), BD3(LDYBD,NZ), BD4(LDXBD,NZ), * BD5(LDXBD,NY), BD6(LDXBD,NY), WORK(NWORK) C C ... LOCAL VARIABLES C LOGICAL SNGULR, PRDX, PRDY, PRDZ INTEGER NXM1, NYM1, NZM1, NML REAL * A, B, C, D, FX, FY, FZ, FYZ, FACTOR, PERTRB, SCALE, UNORM C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 C SNGULR = (BCTY(RIGHT ) .NE. DRCH) .AND. * (BCTY(BOTTOM) .NE. DRCH) .AND. * (BCTY(LEFT ) .NE. DRCH) .AND. * (BCTY(TOP ) .NE. DRCH) .AND. * (BCTY(FRONT ) .NE. DRCH) .AND. * (BCTY(BACK ) .NE. DRCH) .AND. * (ABS(COEFU) .LT. R1MACH(4)) PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC PRDZ = BCTY(FRONT ) .EQ. PRDC C IL = 1 IR = NX JL = 1 JR = NY KL = 1 KR = NZ IF (BCTY(LEFT ) .EQ. DRCH) IL = 2 IF (BCTY(RIGHT ) .EQ. DRCH) IR = NXM1 IF (BCTY(BOTTOM) .EQ. DRCH) JL = 2 IF (BCTY(TOP ) .EQ. DRCH) JR = NYM1 IF (BCTY(BACK ) .EQ. DRCH) KL = 2 IF (BCTY(FRONT ) .EQ. DRCH) KR = NZM1 IF (PRDX) IR = NXM1 IF (PRDY) JR = NYM1 IF (PRDZ) KR = NZM1 NML = (IR-IL+1)*(JR-JL+1)*(KR-KL+1) C C C --------------------------- C CHECK VALIDITY OF ARGUMENTS C --------------------------- C INFO = 0 IF (NX .LT. PRDC) GO TO 901 IF (NY .LT. PRDC) GO TO 902 IF (NZ .LT. PRDC) GO TO 903 IF (LDXU .LT. NX+2) GO TO 904 IF (LDYU .LT. NY+2) GO TO 905 IF ((IORDER .NE. 2) .AND. (IORDER .NE. 4)) GO TO 906 DO 10 K=1,6 IF ((BCTY(K) .LT. DRCH) .OR. (BCTY(K) .GT. PRDC)) GO TO 907 10 CONTINUE IF (((BCTY(RIGHT ) .EQ. PRDC) .AND. (BCTY(LEFT ) .NE. PRDC)) .OR. * ((BCTY(LEFT ) .EQ. PRDC) .AND. (BCTY(RIGHT ) .NE. PRDC)) .OR. * ((BCTY(BOTTOM) .EQ. PRDC) .AND. (BCTY(TOP ) .NE. PRDC)) .OR. * ((BCTY(TOP ) .EQ. PRDC) .AND. (BCTY(BOTTOM) .NE. PRDC)) .OR. * ((BCTY(FRONT ) .EQ. PRDC) .AND. (BCTY(BACK ) .NE. PRDC)) .OR. * ((BCTY(BACK ) .EQ. PRDC) .AND. (BCTY(FRONT ) .NE. PRDC)) ) * GO TO 908 NEEDED = MAX( (NX+1)*(NY+1)*(IORDER-2), * (NX+3)*(NZ+5) + 5*NY + (NX+NZ)/2 + 15 ) IF (NWORK .LT. NEEDED) GO TO 909 IF (LDXGH .LT. NX+1) GO TO 910 IF (LDYGH .LT. NY+1) GO TO 911 C IF (COEFU .GT. 0.0E0) INFO = 1 IF (SNGULR) INFO = 2 C C C ---------------------- C COMPUTE DISCRETIZATION C ---------------------- C IF (IORDER .EQ. 4) THEN CALL HDIS3(NX,NY,NZ,H,COEFU,GH,LDXGH-1,LDYGH-1,BCTY,BD1,BD2, * BD3,BD4,BD5,BD6,LDXBD,LDYBD,A,B,C,D,U,LDXU-1,LDYU-1, * WORK) ELSE CALL FDIS3(NX,NY,NZ,H,COEFU,BCTY,BD1,BD2,BD3,BD4,BD5,BD6, * LDXBD,LDYBD,A,B,C,D,U,LDXU,LDYU) ENDIF C C --------------------------------------- C ADJUST FOR CONSISTENCY IN SINGULAR CASE C --------------------------------------- C PERTRB = 0.0E0 IF (SNGULR) THEN SCALE = 0.0E0 DO 100 K=KL,KR FZ = 1.0E0 IF (.NOT.PRDZ.AND.(K.NE.1).AND.(K.NE.NZ)) FY = 2.0E0 DO 100 J=JL,JR FY = 1.0E0 IF (.NOT.PRDY.AND.(J.NE.1).AND.(J.NE.NY)) FY = 2.0E0 FYZ = FY*FZ DO 100 I=IL,IR FX = 1.0E0 IF (.NOT.PRDX.AND.(I.NE.1).AND.(I.NE.NX)) FX = 2.0E0 FACTOR = FX*FYZ PERTRB = PERTRB + FACTOR*U(I,J,K) SCALE = SCALE + FACTOR 100 CONTINUE PERTRB = -PERTRB/SCALE DO 110 K=KL,KR DO 110 J=JL,JR DO 110 I=IL,IR U(I,J,K) = U(I,J,K) + PERTRB 110 CONTINUE ENDIF C C C ------------------------------ C MATRIX DECOMPOSITION USING FFT C ------------------------------ C CALL MDALG3(A,B,C,D,BCTY,U,LDXU,LDYU,IL,IR,JL,JR,KL,KR,WORK) C C C ----------------------------------------- C SELECT MIN NORM SOLUTION IN SINGULAR CASE C ----------------------------------------- C IF (SNGULR) THEN UNORM = 0.0E0 DO 210 K=KL,KR DO 210 J=JL,JR DO 210 I=IL,IR UNORM = UNORM + U(I,J,K) 210 CONTINUE UNORM = UNORM/REAL(NML) DO 220 K=KL,KR DO 220 J=JL,JR DO 220 I=IL,IR U(I,J,K) = U(I,J,K) - UNORM 220 CONTINUE ENDIF C C C --------------------------------------- C COPY IDENTICAL PLANES IN PERIODIC CASES C --------------------------------------- C IF (PRDX) THEN DO 325 K=1,NZ DO 325 J=1,NY U(NX,J,K) = U(1,J,K) 325 CONTINUE ENDIF C IF (PRDY) THEN DO 375 K=1,NZ DO 375 I=1,NX U(I,NY,K) = U(I,1,K) 375 CONTINUE ENDIF C IF (PRDZ) THEN DO 425 J=1,NY DO 425 I=1,NX U(I,J,NZ) = U(I,J,1) 425 CONTINUE ENDIF C C C ----------- C NORMAL EXIT C ----------- C WORK(1) = PERTRB GO TO 999 C C C ----------- C ERROR EXITS C ----------- C 901 CONTINUE INFO = -1 GO TO 999 C 902 CONTINUE INFO = -2 GO TO 999 C 903 CONTINUE INFO = -3 GO TO 999 C 904 CONTINUE INFO = -4 GO TO 999 C 905 CONTINUE INFO = -5 GO TO 999 C 906 CONTINUE INFO = -6 GO TO 999 C 907 CONTINUE INFO = -7 GO TO 999 C 908 CONTINUE INFO = -8 GO TO 999 C 909 CONTINUE INFO = -9 GO TO 999 C 910 CONTINUE INFO = -10 GO TO 999 C 911 CONTINUE INFO = -11 GO TO 999 C 999 CONTINUE RETURN END SUBROUTINE FDIS3 (NX, NY, NZ, H, COEFU, BCTY, BD1, BD2, BD3, BD4, * BD5, BD6, LDXBD, LDYBD, A, B, C, D, U, * LDXU, LDYU) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C FDIS3 COMPUTES A SECOND ORDER FINITE DIFFERENCE DISCRETIZATION C FOR A THREE-DIMENSIONAL RECTANGULAR DOMAIN C C C P A R A M E T E R S C ------------------- C C NX,NY,NZ INTEGER SCALARS (INPUT) C SEE HFFT3A. C C H, COEFU REAL SCALARS (INPUT) C SEE HFFT3A. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C BD1, BD3 REAL ARRAYS OF SIZE LDYBD BY NZ (INPUT) C SEE HFFT3A. C C BD2, BD4 REAL ARRAYS OF SIZE LDXBD BY NZ (INPUT) C SEE HFFT3A. C C BD5, BD6 REAL ARRAYS OF SIZE LDXBD BY NY (INPUT) C SEE HFFT3A. C C LDXBD INTEGER SCALARS (INPUT) C LDYBD SEE HFFT3A. C C A,B,C,D REAL SCALARS (OUTPUT) C GIVES VALUES IN THE BASIC FINITE DIFFERENCE STENCIL C (SCALED TO O(1)) C C C C C B C C C C C C B C C B A B U = RIGHT SIDE C C B C C C C C C B C C C C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ (INPUT/OUTPUT) C ON INPUT, U(I,J,K) IS THE RIGHT HAND SIDE OF THE PDE C EVALUATED AT THE (I,J,K)TH GRID POINT. C ON OUTPUT, U(I,J,K) IS THE RIGHT HAND SIDE OF THE C DISCRETE PDE AT THE (I,J,K)TH GRID POINT. C C LDXU INTEGER SCALARS (INPUT) C LDYU SEE HFFT3A. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), NX, NY, NZ, LDXU, LDYU REAL * BD1(LDYBD,NZ), BD2(LDXBD,NZ), BD3(LDYBD,NZ), BD4(LDXBD,NZ), * BD5(LDXBD,NY), BD6(LDXBD,NY), U(LDXU,LDYU,*), * COEFU, H, A, B, C, D C C ... LOCAL VARIABLES C LOGICAL PRDX, PRDY, PRDZ, HAVED, HAVEN INTEGER I, J, K, NXM1, NYM1, NZM1 REAL * BETA0, GAMMA0, GEDGE0, GEDGE1, GCORN0, GCORN1 C COMMON /FD3COM/ GAMMA0, GEDGE0, GEDGE1, GCORN0, GCORN1 C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 C HAVED = (BCTY(RIGHT) .EQ. DRCH) .OR. (BCTY(BOTTOM) .EQ. DRCH) .OR. * (BCTY(LEFT ) .EQ. DRCH) .OR. (BCTY(TOP ) .EQ. DRCH) .OR. * (BCTY(FRONT) .EQ. DRCH) .OR. (BCTY(BACK ) .EQ. DRCH) HAVEN = (BCTY(RIGHT) .EQ. NEUM) .OR. (BCTY(BOTTOM) .EQ. NEUM) .OR. * (BCTY(LEFT ) .EQ. NEUM) .OR. (BCTY(TOP ) .EQ. NEUM) .OR. * (BCTY(FRONT) .EQ. NEUM) .OR. (BCTY(BACK ) .EQ. NEUM) PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC PRDZ = BCTY(FRONT ) .EQ. PRDC C H2 = H*H F = -H2*COEFU A = -(24.0E0 + 6.0E0*F) B = 2.0E0 C = 1.0E0 D = 0.0E0 BETA0 = 6.0E0*H2 GAMMA0 = -12.0E0*H GEDGE0 = -10.0E0*H GEDGE1 = -2.0E0*H GCORN0 = -8.0E0*H GCORN1 = -2.0E0*H C C C ---------------------------- C DISCRETIZE RIGHT SIDE OF PDE C ---------------------------- C ISTRT = 2 ISTOP = NXM1 IF (BCTY(LEFT ) .NE. DRCH) ISTRT = 1 IF (BCTY(RIGHT ) .EQ. NEUM) ISTOP = NX JSTRT = 2 JSTOP = NYM1 IF (BCTY(BOTTOM) .NE. DRCH) JSTRT = 1 IF (BCTY(TOP ) .EQ. NEUM) JSTOP = NY KSTRT = 2 KSTOP = NZM1 IF (BCTY(BACK ) .NE. DRCH) KSTRT = 1 IF (BCTY(FRONT ) .EQ. NEUM) KSTOP = NZ C DO 100 K=KSTRT,KSTOP DO 100 J=JSTRT,JSTOP DO 100 I=ISTRT,ISTOP U(I,J,K) = BETA0*U(I,J,K) 100 CONTINUE C C C --------------------------------------- C DISCRETIZE POINTS ON NEUMANN BOUNDARIES C --------------------------------------- C IF (HAVEN) CALL FD3N(NX,NY,NZ,BCTY,BD1,BD2,BD3,BD4,BD5,BD6, * LDXBD,LDYBD,U,LDXU,LDYU) C C ---------------------------------------------- C ADJUST POINTS ADJACENT TO DIRICHLET BOUNDARIES C ---------------------------------------------- C IF (HAVED) THEN CALL STORD3(NX,NY,NZ,BCTY,BD1,BD2,BD3,BD4,BD5,BD6,LDXBD,LDYBD, * U,LDXU,LDYU) CALL FD3D(NX,NY,NZ,BCTY,B,C,D,U,LDXU,LDYU) ENDIF C C C ---- C EXIT C ---- C RETURN END SUBROUTINE FD2DA (NX, NY, UD, LDXU, BCTY, A, B, C, U) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C FD2DA ELIMINATES TERMS FROM EQUATIONS CORRESPONDING TO ALL POINTS C OF A TWO-DIMENSIONAL RECTANGULAR DOMAIN ACCORDING TO A GIVEN C STENCIL. C C (THIS ROUTINE IS A UTILITY USED TO ELIMINATE AN ENTIRE DIRICHLET C PLANE IN THREE-DIMENSIONAL CALCULATIONS.) C C C P A R A M E T E R S C ------------------- C C NX, NY INTEGER SCALARS (INPUT) C SEE HFFT2A. C C UD REAL ARRAY OF SIZE LDXU BY NY (INPUT) C ALL ENTRIES CONTAIN VALUES OF THE FUNCTION TO BE C ELIMINATED. C C LDXU INTEGER SCALAR (INPUT) C THE LEADING DIMENSION OF THE ARRAYS UD AND U EXACTLY C AS DECLARED IN THE CALLING PROGRAM. C C BCTY INTEGER ARRAY OF SIZE 4 (INPUT) C SEE HFFT2A. C C A,B,C REAL SCALARS (INPUT) C FINITE DIFFERENCE STENCIL COEFFICIENTS TO USE IN THE C ELIMINATION C C C B C C B A B C C B C C C U REAL ARRAY OF SIZE LDXU BY NY (INPUT/OUTPUT) C ON INPUT, ENTRIES CORRESPONDING TO POINTS WHERE THE C SOLUTION IS TO BE DETERMINED CONTAIN THE RIGHT HAND C SIDE OF A FINITE DIFFERENCE DISCRETIZATION. C ON EXIT, THESE ENTRIES ARE UPDATED SUCH THAT C THE GRID FUNCTION IN UD IS SUBTRACTED USING THE C GIVEN STENCIL. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(4), NX, NY REAL * UD(LDXU,NY), U(LDXU,NY), A, B, C C C ... LOCAL VARIABLES C LOGICAL PRDX, PRDY INTEGER I, J, NXM1, NYM1 REAL * TWOB, TWOC C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NXM1 = NX - 1 NYM1 = NY - 1 C PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC C TWOB = 2.0E0*B TWOC = 2.0E0*C FORC = 4.0E0*C C C C ----------------------- C PROCESS INTERIOR POINTS C ----------------------- C DO 100 J=2,NYM1 DO 100 I=2,NXM1 U(I,J) = U(I,J) - (C*UD(I-1,J+1)+B*UD(I,J+1)+C*UD(I+1,J+1) * +B*UD(I-1,J )+A*UD(I,J )+B*UD(I+1,J ) * +C*UD(I-1,J-1)+B*UD(I,J-1)+C*UD(I+1,J-1)) 100 CONTINUE C C C ---------------------- C PROCESS NEUMANN POINTS C ---------------------- C C ... NEUMANN POINTS ON RIGHT EDGE C IF (BCTY(RIGHT) .EQ. NEUM) THEN DO 150 J=2,NYM1 U(NX,J) = U(NX,J) - (TWOC*UD(NXM1,J+1)+B*UD(NX,J+1) * +TWOB*UD(NXM1,J )+A*UD(NX,J ) * +TWOC*UD(NXM1,J-1)+B*UD(NX,J-1)) 150 CONTINUE ENDIF C C ... NEUMANN POINTS ON LEFT EDGE C IF (BCTY(LEFT) .EQ. NEUM) THEN DO 250 J=2,NYM1 U(1,J) = U(1,J) - (B*UD(1,J+1)+TWOC*UD(2,J+1) * +A*UD(1,J )+TWOB*UD(2,J ) * +B*UD(1,J-1)+TWOC*UD(2,J-1)) 250 CONTINUE ENDIF C C ... NEUMANN POINTS ON BOTTOM EDGE C IF (BCTY(BOTTOM) .EQ. NEUM) THEN DO 350 I=2,NXM1 U(I,1) = U(I,1) - (TWOC*UD(I+1,2)+B*UD(I+1,1) * +TWOB*UD(I ,2)+A*UD(I ,1) * +TWOC*UD(I-1,2)+B*UD(I-1,1)) 350 CONTINUE ENDIF C C ... NEUMANN POINTS ON TOP EDGE C IF (BCTY(TOP) .EQ. NEUM) THEN DO 450 I=2,NXM1 U(I,NY) = U(I,NY) - (TWOC*UD(I+1,NYM1)+B*UD(I+1,NY) * +TWOB*UD(I ,NYM1)+A*UD(I ,NY) * +TWOC*UD(I-1,NYM1)+B*UD(I-1,NY)) 450 CONTINUE ENDIF C C ... BOTTOM RIGHT NEUMANN CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) * U(NX,1) = U(NX,1) - (TWOB*UD(NX ,2)+ A*UD(NX ,1) * +FORC*UD(NXM1,2)+TWOB*UD(NXM1,1)) C C ... BOTTOM LEFT NEUMANN CORNER C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) * U(1,1) = U(1,1) - (TWOB*UD(1,2)+ A*UD(1,1) * +FORC*UD(2,2)+TWOB*UD(2,1)) C C ... TOP RIGHT NEUMANN CORNER C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(RIGHT) .EQ. NEUM)) * U(NX,NY) = U(NX,NY) - (TWOB*UD(NXM1,NY )+ A*UD(NX,NY ) * +FORC*UD(NXM1,NYM1)+TWOB*UD(NX,NYM1)) C C ... TOP LEFT NEUMANN CORNER C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(LEFT) .EQ. NEUM)) * U(1,NY) = U(1,NY) - (TWOB*UD(2,NY )+ A*UD(1,NY ) * +FORC*UD(2,NYM1)+TWOB*UD(1,NYM1)) C C C ----------------------- C PROCESS PERIODIC POINTS C ----------------------- C C ... LEFT PERIODIC EDGE C IF (PRDX) THEN DO 605 J=2,NYM1 U(1,J) = U(1,J) - (C*UD(NXM1,J+1)+B*UD(1,J+1)+C*UD(2,J+1) * +B*UD(NXM1,J )+A*UD(1,J )+B*UD(2,J ) * +C*UD(NXM1,J-1)+B*UD(1,J-1)+C*UD(2,J-1)) 605 CONTINUE ENDIF C C ... BOTTOM PERIODIC EDGE C IF (PRDY) THEN DO 615 I=2,NXM1 U(I,1) = U(I,1)-(C*UD(I-1,2 )+B*UD(I,2 )+C*UD(I+1,2 ) * +B*UD(I-1,1 )+A*UD(I,1 )+B*UD(I+1,1 ) * +C*UD(I-1,NYM1)+B*UD(I,NYM1)+C*UD(I+1,NYM1)) 615 CONTINUE ENDIF C C ... LEFT BOTTOM PERIODIC CORNER C IF (PRDX .AND. PRDY) * U(1,1) = U(1,1) - (C*UD(NXM1,2 )+B*UD(1,2 )+C*UD(2,2 ) * +B*UD(NXM1,1 )+A*UD(1,1 )+B*UD(2,1 ) * +C*UD(NXM1,NYM1)+B*UD(1,NYM1)+C*UD(2,NYM1)) C C ... LEFT BOTTOM PERIODIC/NEUMANN CORNER C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM)) * U(1,1) = U(1,1) - (TWOC*UD(2 ,2)+B*UD(2 ,1) * +TWOB*UD(1 ,2)+A*UD(1 ,1) * +TWOC*UD(NXM1,2)+B*UD(NXM1,1)) C C ... LEFT BOTTOM NEUMANN/PERIODIC CORNER C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY) * U(1,1) = U(1,1) - (B*UD(1,2 )+TWOC*UD(2,2 ) * +A*UD(1,1 )+TWOB*UD(2,1 ) * +B*UD(1,NYM1)+TWOC*UD(2,NYM1)) C C ... LEFT TOP PERIODIC/NEUMANN CORNER C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM)) * U(1,NY) = U(1,NY) - (TWOC*UD(2 ,NYM1)+B*UD(2 ,NY) * +TWOB*UD(1 ,NYM1)+A*UD(1 ,NY) * +TWOC*UD(NXM1,NYM1)+B*UD(NXM1,NY)) C C ... RIGHT BOTTOM NEUMANN/PERIODIC CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY) * U(NX,1) = U(NX,1) - (TWOC*UD(NXM1,2 )+B*UD(NX,2 ) * +TWOB*UD(NXM1,1 )+A*UD(NX,1 ) * +TWOC*UD(NXM1,NYM1)+B*UD(NX,NYM1)) C C C ---- C EXIT C ---- C RETURN END SUBROUTINE FD3D (NX, NY, NZ, BCTY, B, C, D, U, LDXU, LDYU) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C FD3D ELIMINATES KNOWN TERMS FROM EQUATIONS CORRESPONDING TO POINTS C NEAR DIRICHLET BOUNDARIES OF A THREE-DIMENSIONAL RECTANGULAR DOMAIN. C C C P A R A M E T E R S C ------------------- C C NX,NY,NZ INTEGER SCALARS (INPUT) C SEE HFFT3A. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C UD REAL ARRAY OF SIZE LDXU BY LDYU BY NZ (INPUT) C ENTRIES CORRESPONDING TO DIRICHLET POINTS CONTAIN C KNOWN VALUES OF THE SOLUTION. C C B,C,D REAL SCALARS (INPUT) C FINITE DIFFERENCE STENCIL COEFFICIENTS TO USE IN THE C ELIMINATION C C D C D C C B C C D C D C C C B C C B * B C C B C C C D C D C C B C C D C D C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ (INPUT/OUTPUT) C ON INPUT, ENTRIES CORRESPONDING TO POINTS WHERE THE C SOLUTION IS TO BE DETERMINED CONTAIN THE RIGHT HAND C SIDE OF A FINITE DIFFERENCE DISCRETIZATION. C ON EXIT, THESE ENTRIES ARE UPDATED SUCH THAT KNOWN C TERMS (DIRICHLET POINTS) ARE ELIMINATED FROM THE C LEFT HAND SIDE OF THE EQUATION. C C LDXU INTEGER SCALAR (INPUT) C THE LEADING DIMENSION OF THE ARRAYS UD AND U EXACTLY C AS DECLARED IN THE CALLING PROGRAM. C C LDYU INTEGER SCALAR (INPUT) C THE SECOND DIMENSION OF THE ARRAYS UD AND U EXACTLY C AS DECLARED IN THE CALLING PROGRAM. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), NX, NY, NZ REAL * U(LDXU,LDYU,NZ), B, C, D C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NZM1 = NZ - 1 NZM2 = NZ - 2 C TWOC = 2.0E0*C TWOD = 2.0E0*D C C ------------------------ C PROCESS EACH (X,Y) PLANE C ------------------------ C C ... BACK PLANE C IF (BCTY(BACK) .EQ. NEUM) THEN CALL FD2D(NX,NY,U(1,1,1),LDXU,BCTY,B,C,U(1,1,1)) CALL FD2D(NX,NY,U(1,1,2),LDXU,BCTY,TWOC,TWOD,U(1,1,1)) ENDIF IF (BCTY(BACK) .EQ. PRDC) THEN CALL FD2D(NX,NY,U(1,1,NZM1),LDXU,BCTY,C,D,U(1,1,1)) CALL FD2D(NX,NY,U(1,1,1),LDXU,BCTY,B,C,U(1,1,1)) CALL FD2D(NX,NY,U(1,1,2),LDXU,BCTY,C,D,U(1,1,1)) ENDIF C C ... INTERIOR PLANE NEAR BACK C IF (BCTY(BACK) .EQ. DRCH) * CALL FD2DA(NX,NY,U(1,1,1),LDXU,BCTY,B,C,D,U(1,1,2)) IF (BCTY(BACK) .NE. DRCH) * CALL FD2D(NX,NY,U(1,1,1),LDXU,BCTY,C,D,U(1,1,2)) CALL FD2D(NX,NY,U(1,1,2),LDXU,BCTY,B,C,U(1,1,2)) CALL FD2D(NX,NY,U(1,1,3),LDXU,BCTY,C,D,U(1,1,2)) C C ... INTERIOR PLANES C DO 100 K=3,NZM2 CALL FD2D(NX,NY,U(1,1,K-1),LDXU,BCTY,C,D,U(1,1,K)) CALL FD2D(NX,NY,U(1,1,K ),LDXU,BCTY,B,C,U(1,1,K)) CALL FD2D(NX,NY,U(1,1,K+1),LDXU,BCTY,C,D,U(1,1,K)) 100 CONTINUE C C ... INTERIOR PLANE NEAR FRONT C CALL FD2D(NX,NY,U(1,1,NZM2),LDXU,BCTY,C,D,U(1,1,NZM1)) CALL FD2D(NX,NY,U(1,1,NZM1),LDXU,BCTY,B,C,U(1,1,NZM1)) IF (BCTY(FRONT) .EQ. DRCH) * CALL FD2DA(NX,NY,U(1,1,NZ),LDXU,BCTY,B,C,D,U(1,1,NZM1)) IF (BCTY(FRONT) .EQ. NEUM) * CALL FD2D(NX,NY,U(1,1,NZ),LDXU,BCTY,C,D,U(1,1,NZM1)) IF (BCTY(FRONT) .EQ. PRDC) * CALL FD2D(NX,NY,U(1,1,1),LDXU,BCTY,C,D,U(1,1,NZM1)) C C ... FRONT PLANE C IF (BCTY(FRONT) .EQ. NEUM) THEN CALL FD2D(NX,NY,U(1,1,NZM1),LDXU,BCTY,TWOC,TWOD,U(1,1,NZ)) CALL FD2D(NX,NY,U(1,1,NZ ),LDXU,BCTY,B,C,U(1,1,NZ)) ENDIF C C C ---- C EXIT C ---- C RETURN END SUBROUTINE FD3N (NX, NY, NZ, BCTY, BD1, BD2, BD3, BD4, BD5, BD6, * LDXBD, LDYBD, U, LDXU, LDYU) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C FD3N COMPUTES THE SECOND ORDER FINITE DIFFERENCE DISCRETIZATION C AT ALL BOUNDARY POINTS OF A THREE-DIMENSIONAL RECTANGULAR DOMAIN C WHERE NEUMANN BOUNDARY CONDITIONS HAVE BEEN SPECIFIED. C C C P A R A M E T E R S C ------------------- C C NX,NY,NZ INTEGER SCALARS (INPUT) C SEE HFFT3A. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C BD1, BD3 REAL ARRAYS OF SIZE LDYBD BY NZ (INPUT) C SEE HFFT3A. C C BD2, BD4 REAL ARRAYS OF SIZE LDXBD BY NZ (INPUT) C SEE HFFT3A. C C BD5, BD6 REAL ARRAYS OF SIZE LDXBD BY NY (INPUT) C SEE HFFT3A. C C LDXBD INTEGER SCALARS (INPUT) C LDYBD SEE HFFT3A. C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ (OUTPUT) C ON EXIT, ENTRIES OF U CORRESPONDING TO NEUMANN BOUNDARY C POINTS CONTAIN THE RIGHT HAND SIDE OF THE SECOND ORDER C FINITE DIFFERENCE DISCRETIZATION C C LDXU INTEGER SCALARS (INPUT) C LDYU SEE HFFT3A. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), NX, NY, NZ, LDXBD, LDYBD, LDXU, LDYU REAL * U(LDXU,LDYU,NZ), BD1(LDYBD,NZ), BD2(LDXBD,NZ), BD3(LDYBD,NZ), * BD4(LDXBD,NZ), BD5(LDXBD,NY), BD6(LDXBD,NY) C C ... LOCAL VARIABLES C LOGICAL PRDX, PRDY, PRDZ INTEGER I, J, K, NXM1, NYM1, NZM1 REAL * GAMMA0, GEDGE0, GEDGE1, GCORN0, GCORN1 C COMMON /FD3COM/ GAMMA0, GEDGE0, GEDGE1, GCORN0, GCORN1 C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC PRDZ = BCTY(FRONT ) .EQ. PRDC C C C ------------------------------ C HANDLE POINTS ON NEUMANN SIDES C ------------------------------ C C ... RIGHT SIDE C IF (BCTY(RIGHT) .EQ. NEUM) THEN DO 205 K=2,NZM1 DO 205 J=2,NYM1 U(NX,J,K) = U(NX,J,K) + GAMMA0*BD1(J,K) 205 CONTINUE ENDIF C C ... LEFT SIDE C IF (BCTY(LEFT) .EQ. NEUM) THEN DO 215 K=2,NZM1 DO 215 J=2,NYM1 U(1,J,K) = U(1,J,K) - GAMMA0*BD3(J,K) 215 CONTINUE ENDIF C C ... TOP SIDE C IF (BCTY(TOP) .EQ. NEUM) THEN DO 225 K=2,NZM1 DO 225 I=2,NXM1 U(I,NY,K) = U(I,NY,K) + GAMMA0*BD4(I,K) 225 CONTINUE ENDIF C C ... BOTTOM SIDE C IF (BCTY(BOTTOM) .EQ. NEUM) THEN DO 235 K=2,NZM1 DO 235 I=2,NXM1 U(I,1,K) = U(I,1,K) - GAMMA0*BD2(I,K) 235 CONTINUE ENDIF C C ... FRONT SIDE C IF (BCTY(FRONT) .EQ. NEUM) THEN DO 245 J=2,NYM1 DO 245 I=2,NXM1 U(I,J,NZ) = U(I,J,NZ) + GAMMA0*BD5(I,J) 245 CONTINUE ENDIF C C ... BACK SIDE C IF (BCTY(BACK) .EQ. NEUM) THEN DO 255 J=2,NYM1 DO 255 I=2,NXM1 U(I,J,1) = U(I,J,1) - GAMMA0*BD6(I,J) 255 CONTINUE ENDIF C C C ------------------------------ C HANDLE POINTS ON NEUMANN EDGES C ------------------------------ C C ... RIGHT TOP EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 265 K=2,NZM1 U(NX,NY,K) = U(NX,NY,K) * + GEDGE0*( BD1(NY,K) + BD4(NX,K) ) * + GEDGE1*( BD1(NYM1,K) + BD4(NXM1,K) ) 265 CONTINUE ENDIF C C ... RIGHT BOTTOM EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) THEN DO 275 K=2,NZM1 U(NX,1,K) = U(NX,1,K) * + GEDGE0*( BD1(1,K) - BD2(NX,K) ) * + GEDGE1*( BD1(2,K) - BD2(NXM1,K) ) 275 CONTINUE ENDIF C C ... LEFT TOP EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 285 K=2,NZM1 U(1,NY,K) = U(1,NY,K) * + GEDGE0*( -BD3(NY,K) + BD4(1,K) ) * + GEDGE1*( -BD3(NYM1,K) + BD4(2,K) ) 285 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) THEN DO 295 K=2,NZM1 U(1,1,K) = U(1,1,K) * + GEDGE0*( -BD3(1,K) - BD2(1,K) ) * + GEDGE1*( -BD3(2,K) - BD2(2,K) ) 295 CONTINUE ENDIF C C ... RIGHT FRONT EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 305 J=2,NYM1 U(NX,J,NZ) = U(NX,J,NZ) * + GEDGE0*( BD1(J,NZ) + BD5(NX,J) ) * + GEDGE1*( BD1(J,NZM1) + BD5(NXM1,J) ) 305 CONTINUE ENDIF C C ... RIGHT BACK EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 315 J=2,NYM1 U(NX,J,1) = U(NX,J,1) * + GEDGE0*( BD1(J,1) - BD6(NX,J) ) * + GEDGE1*( BD1(J,2) - BD6(NXM1,J) ) 315 CONTINUE ENDIF C C ... LEFT FRONT EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 325 J=2,NYM1 U(1,J,NZ) = U(1,J,NZ) * + GEDGE0*( -BD3(J,NZ) + BD5(1,J) ) * + GEDGE1*( -BD3(J,NZM1) + BD5(2,J) ) 325 CONTINUE ENDIF C C ... LEFT BACK EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 335 J=2,NYM1 U(1,J,1) = U(1,J,1) * + GEDGE0*( -BD3(J,1) - BD6(1,J) ) * + GEDGE1*( -BD3(J,2) - BD6(2,J) ) 335 CONTINUE ENDIF C C ... TOP FRONT EDGE C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 345 I=2,NXM1 U(I,NY,NZ) = U(I,NY,NZ) * + GEDGE0*( BD4(I,NZ) + BD5(I,NY) ) * + GEDGE1*( BD4(I,NZM1) + BD5(I,NYM1) ) 345 CONTINUE ENDIF C C ... TOP BACK EDGE C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 355 I=2,NXM1 U(I,NY,1) = U(I,NY,1) * + GEDGE0*( BD4(I,1) - BD6(I,NY) ) * + GEDGE1*( BD4(I,2) - BD6(I,NYM1) ) 355 CONTINUE ENDIF C C ... BOTTOM FRONT EDGE C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 365 I=2,NXM1 U(I,1,NZ) = U(I,1,NZ) * + GEDGE0*( -BD2(I,NZ) + BD5(I,1) ) * + GEDGE1*( -BD2(I,NZM1) + BD5(I,2) ) 365 CONTINUE ENDIF C C ... BOTTOM BACK EDGE C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 375 I=2,NXM1 U(I,1,1) = U(I,1,1) * + GEDGE0*( -BD2(I,1) - BD6(I,1) ) * + GEDGE1*( -BD2(I,2) - BD6(I,2) ) 375 CONTINUE ENDIF C C C -------------------------------- C HANDLE POINTS AT NEUMANN CORNERS C -------------------------------- C C ... RIGHT TOP FRONT CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(NX,NY,NZ) = U(NX,NY,NZ) * + GCORN0*( BD1(NY,NZ)+BD4(NX,NZ)+BD5(NX,NY) ) * + GCORN1*( BD1(NYM1,NZ) + BD1(NY,NZM1) * + BD4(NXM1,NZ) + BD4(NX,NZM1) * + BD5(NXM1,NY) + BD5(NX,NYM1) ) C C ... LEFT TOP FRONT CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(1,NY,NZ) = U(1,NY,NZ) * + GCORN0*(-BD3(NY,NZ)+BD4(1,NZ)+BD5(1,NY) ) * + GCORN1*(-BD3(NYM1,NZ) - BD3(NY,NZM1) * + BD4(2,NZ) + BD4(1,NZM1) * + BD5(2,NY) + BD5(1,NYM1) ) C C ... RIGHT BOTTOM FRONT CORNER C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM)) * U(NX,1,NZ) = U(NX,1,NZ) * + GCORN0*( BD1(1,NZ)-BD2(NX,NZ)+BD5(NZ,1) ) * + GCORN1*( BD1(2,NZ) + BD1(1,NZM1) * - BD2(NXM1,NZ) - BD2(NX,NZM1) * + BD5(NXM1,1) + BD5(NX,2) ) C C ... LEFT BOTTOM FRONT CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM)) * U(1,1,NZ) = U(1,1,NZ) * + GCORN0*(-BD3(1,NZ)-BD2(1,NZ)+BD5(1,1) ) * + GCORN1*(-BD3(2,NZ) - BD3(1,NZM1) * - BD2(2,NZ) - BD2(1,NZM1) * + BD5(2,1) + BD5(1,2) ) C C ... RIGHT TOP BACK CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(NX,NY,1) = U(NX,NY,1) * + GCORN0*( BD1(NY,1)+BD4(NX,1)-BD6(NX,NY) ) * + GCORN1*( BD1(NYM1,1) + BD1(NY,2) * + BD4(NXM1,1) + BD4(NX,2) * - BD6(NXM1,NY) - BD6(NX,NYM1) ) C C ... LEFT TOP BACK CORNER C IF ((BCTY(LEFT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(BACK) .EQ. NEUM)) * U(1,NY,1) = U(1,NY,1) * + GCORN0*(-BD3(NY,1)+BD4(1,1)-BD6(1,NY) ) * + GCORN1*(-BD3(NYM1,1) - BD3(NY,2) * + BD4(2,1) + BD4(1,2) * - BD6(2,NY) - BD6(1,NYM1) ) C C ... RIGHT BOTTOM BACK CORNER C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(NX,1,1) = U(NX,1,1) * + GCORN0*( BD1(1,1)-BD2(NX,1)-BD6(NX,1) ) * + GCORN1*( BD1(2,1) + BD1(1,2) * - BD2(NXM1,1) - BD2(NX,2) * - BD6(NXM1,1) - BD6(NX,2) ) C C ... LEFT BOTTOM BACK CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(1,1,1) = U(1,1,1) * + GCORN0*(-BD3(1,1)-BD2(1,1)-BD6(1,1) ) * + GCORN1*(-BD3(2,1) - BD3(1,2) * - BD2(2,1) - BD2(1,2) * - BD6(2,1) - BD6(1,2) ) C C IF (.NOT. (PRDX .OR. PRDY .OR. PRDZ)) GO TO 999 C C C --------------------------------------- C HANDLE POINTS AT NEUMANN/PERIODIC EDGES C --------------------------------------- C C ... BOTTOM BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) THEN DO 405 I=2,NXM1 U(I,1,1) = U(I,1,1) - GAMMA0*BD2(I,1) 405 CONTINUE ENDIF C C ... BOTTOM BACK EDGE (PERIODIC/NEUMANN) C IF (PRDY .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 415 I=2,NXM1 U(I,1,1) = U(I,1,1) - GAMMA0*BD6(I,1) 415 CONTINUE ENDIF C C ... LEFT BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDZ) THEN DO 425 J=2,NYM1 U(1,J,1) = U(1,J,1) - GAMMA0*BD3(J,1) 425 CONTINUE ENDIF C C ... LEFT BACK EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 435 J=2,NYM1 U(1,J,1) = U(1,J,1) - GAMMA0*BD6(1,J) 435 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE (PERIODIC/NEUNANN) C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM)) THEN DO 445 K=2,NZM1 U(1,1,K) = U(1,1,K) - GAMMA0*BD2(1,K) 445 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE (NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY) THEN DO 455 K=2,NZM1 U(1,1,K) = U(1,1,K) - GAMMA0*BD3(1,K) 455 CONTINUE ENDIF C C ... LEFT TOP EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 465 K=2,NZM1 U(1,NY,K) = U(1,NY,K) + GAMMA0*BD4(1,K) 465 CONTINUE ENDIF C C ... LEFT FRONT EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 475 J=2,NYM1 U(1,J,NZ) = U(1,J,NZ) + GAMMA0*BD5(1,J) 475 CONTINUE ENDIF C C ... RIGHT BOTTOM EDGE (NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY) THEN DO 485 K=2,NZM1 U(NX,1,K) = U(NX,1,K) + GAMMA0*BD1(1,K) 485 CONTINUE ENDIF C C ... BOTTOM FRONT EDGE (PERIODIC/NEUMANN) C IF (PRDY .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 495 I=2,NXM1 U(I,1,NZ) = U(I,1,NZ) + GAMMA0*BD5(I,1) 495 CONTINUE ENDIF C C ... RIGHT BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDZ) THEN DO 505 J=2,NYM1 U(NX,J,1) = U(NX,J,1) + GAMMA0*BD1(J,1) 505 CONTINUE ENDIF C C ... TOP BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(TOP) .EQ. NEUM) .AND. PRDZ) THEN DO 515 I=2,NXM1 U(I,NY,1) = U(I,NY,1) + GAMMA0*BD4(I,1) 515 CONTINUE ENDIF C C C ----------------------------------------- C HANDLE POINTS AT NEUMANN/PERIODIC CORNERS C ----------------------------------------- C C ... LEFT BOTTOM BACK CORNER (NEUMANN/PERIODIC/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY .AND. PRDZ) * U(1,1,1) = U(1,1,1) - GAMMA0*BD3(1,1) C C ... LEFT BOTTOM BACK CORNER (PERIODIC/NEUMANN/PERIODIC) C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) * U(1,1,1) = U(1,1,1) - GAMMA0*BD2(1,1) C C ... LEFT BOTTOM BACK CORNER (PERIODIC/PERIODIC/NEUMANN) C IF (PRDX .AND. PRDY .AND. (BCTY(BACK) .EQ. NEUM)) * U(1,1,1) = U(1,1,1) - GAMMA0*BD6(1,1) C C ... LEFT BOTTOM BACK CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD2(1,1) - BD6(1,1) ) * + GEDGE1*( -BD2(1,2) - BD6(1,2) ) C C ... LEFT BOTTOM BACK CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(BACK) .EQ. NEUM)) * U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD3(1,1) - BD6(1,1) ) * + GEDGE1*( -BD3(1,2) - BD6(2,1) ) C C ... LEFT BOTTOM BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) * U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD3(1,1) - BD2(1,1) ) * + GEDGE1*( -BD3(2,1) - BD2(2,1) ) C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/PERIODIC/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. PRDZ) * U(NX,1,1) = U(NX,1,1) + GAMMA0*BD1(1,1) C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) * U(NX,1,1) = U(NX,1,1) * + GEDGE0*( BD1(1,1) - BD2(NX,1) ) * + GEDGE1*( BD1(2,1) - BD2(NXM1,1) ) C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(NX,1,1) = U(NX,1,1) * + GEDGE0*( BD1(1,1) - BD6(NX,1) ) * + GEDGE1*( BD1(1,2) - BD6(NXM1,1) ) C C ... LEFT TOP BACK CORNER (PERIODIC/NEUMANN/PERIODIC) C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM) .AND. PRDZ) * U(1,NY,1) = U(1,NY,1) + GAMMA0*BD4(1,1) C C ... LEFT TOP BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. PRDZ) * U(1,NY,1) = U(1,NY,1) * + GEDGE0*( -BD3(NY,1) + BD4(1,1) ) * + GEDGE1*( -BD3(NYM1,1) + BD4(2,1) ) C C ... LEFT TOP BACK CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(BACK) .EQ. NEUM)) * U(1,NY,1) = U(1,NY,1) * + GEDGE0*( BD4(1,1) - BD6(1,NY) ) * + GEDGE1*( BD4(1,2) - BD6(1,NYM1) ) C C ... LEFT BOTTOM FRONT CORNER (PERIODIC/PERIODIC/NEUMANN) C IF (PRDX .AND. PRDY .AND. (BCTY(FRONT) .EQ. NEUM)) * U(1,1,NZ) = U(1,1,NZ) + GAMMA0*BD5(1,1) C C ... LEFT BOTTOM FRONT CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(1,1,NZ) = U(1,1,NZ) * + GEDGE0*( -BD3(1,NZ) + BD5(1,1) ) * + GEDGE1*( -BD3(1,NZM1) + BD5(2,1) ) C C ... LEFT BOTTOM FRONT CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM)) * U(1,1,NZ) = U(1,1,NZ) * + GEDGE0*( -BD2(1,NZ) + BD5(1,1) ) * + GEDGE1*( -BD2(1,NZM1) + BD5(1,2) ) C C ... LEFT TOP FRONT CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(1,NY,NZ) = U(1,NY,NZ) * + GEDGE0*( BD4(1,NZ) + BD5(1,NY) ) * + GEDGE1*( BD4(1,NZM1) + BD5(1,NYM1) ) C C ... RIGHT BOTTOM FRONT CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(NX,1,NZ) = U(NX,1,NZ) * + GEDGE0*( BD1(1,NZ) + BD5(NX,1) ) * + GEDGE1*( BD1(1,NZM1) + BD5(NXM1,1) ) C C ... RIGHT TOP BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. PRDZ) * U(NX,NY,1) = U(NX,NY,1) * + GEDGE0*( BD1(NY,1) + BD4(NX,1) ) * + GEDGE1*( BD1(NYM1,1) + BD4(NXM1,1) ) C C C ---- C EXIT C ---- C 999 CONTINUE RETURN END SUBROUTINE HDIS3 (NX, NY, NZ, H, COEFU, GH, LMXGH, LMYGH, BCTY, * BD1, BD2, BD3, BD4, BD5, BD6, LDXBD, LDYBD, * A, B, C, D, U, LMXU, LMYU, WPLANE) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C HDIS3 COMPUTES A FOURTH ORDER FINITE DIFFERENCE DISCRETIZATION C FOR A THREE-DIMENSIONAL RECTANGULAR DOMAIN C C C P A R A M E T E R S C ------------------- C C NX,NY,NZ INTEGER SCALARS (INPUT) C SEE HFFT3A. C C H, COEFU REAL SCALARS (INPUT) C SEE HFFT3A. C C GH REAL ARRAY OF SIZE LMXGH+1 BY LMYGH+1 BY NZ+1 (INPUT) C SEE HFFT3A. C C LMXGH INTEGER SCALAR (INPUT) C UPPER LIMIT OF FIRST DIMENSION OF ARRAY GH. MUST BE C LDXGH-1 (LDXGH IS DEFINED IN HFFT3A). C C LMYGH INTEGER SCALAR (INPUT) C UPPER LIMIT OF SECOND DIMENSION OF ARRAY GH. MUST BE C LDYGH-1 (LDYGH IS DEFINED IN HFFT3A). C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C BD1, BD3 REAL ARRAYS OF SIZE LDYBD BY NZ (INPUT) C SEE HFFT3A. C C BD2, BD4 REAL ARRAYS OF SIZE LDXBD BY NZ (INPUT) C SEE HFFT3A. C C BD5, BD6 REAL ARRAYS OF SIZE LDXBD BY NY (INPUT) C SEE HFFT3A. C C LDXBD INTEGER SCALARS (INPUT) C LDYBD SEE HFFT3A. C C C A,B,C,D REAL SCALARS (OUTPUT) C GIVES VALUES IN THE BASIC FINITE DIFFERENCE STENCIL C (SCALED TO O(1)) C C C C C B C C C C C C B C C B A B U = RIGHT SIDE C C B C C C C C C B C C C C C U REAL ARRAY OF SIZE LMXU+1 BY LMYU+1 BY NZ+2 C (INPUT/OUTPUT) C ON INPUT, U(I,J,K) IS THE RIGHT HAND SIDE OF THE PDE C EVALUATED AT THE (I,J,K)TH GRID POINT. C ON OUTPUT, U(I,J,K) IS THE RIGHT HAND SIDE OF THE C DISCRETE PDE AT THE (I,J,K)TH GRID POINT. C C LMXU INTEGER SCALAR (INPUT) C UPPER LIMIT OF FIRST DIMENSION OF ARRAY U. MUST BE C LDXU-1 (LDXU IS DEFINED IN HFFT3A). C C LMYU INTEGER SCALAR (INPUT) C UPPER LIMIT OF SECOND DIMENSION OF ARRAY U. MUST BE C LDYU-1 (LDYU IS DEFINED IN HFFT3A). C C WPLANE REAL ARRAY OF SIZE NX+1 BY 2 C WORKING STORAGE FOR HDIS3. C C C ****************************************************************** C * * C * NOTE -- THE ARRAYS U AND GH ARE INDEXED DIFFERENTLY IN * C * THIS ROUTINE: U(0:LMXU,0:LMYU,0:*) AND * C * GH(0:LMXGH,0:LMYGH,0:*) * C * * C ****************************************************************** C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), NX, NY, NZ, LMXGH, LMYGH, LDXBD, LDYBD, * LMXU, LMYU REAL * GH(0:LMXGH,0:LMYGH,0:*), U(0:LMXU,0:LMYU,0:*), BD1(LDYBD,NZ), * BD2(LDXBD,NZ), BD3(LDYBD,NZ), BD4(LDXBD,NZ), BD5(LDXBD,NY), * BD6(LDXBD,NY), WPLANE(0:NX,0:NY,0:1), * COEFU, H, A, B, C, D C C ... LOCAL VARIABLES C LOGICAL PRDX, PRDY, PRDZ, HELMHZ, HAVED, HAVEN INTEGER I, J, K, NXM1, NYM1, NZM1, LDXU, LDYU, P0, PM1 REAL * F, F2, H2, BETA0, BETA1, BETA2, GAMMA0, GAMMA1, * GEDGE0, GEDGE1, GEDGE2, GEDGE3, GEDGE4, GEDGE5, * GCORN0, GCORN1, GCORN2, GCORN3, GCORN4, GCORN5 C COMMON /HD3COM/ GAMMA0, GAMMA1, * GEDGE0, GEDGE1, GEDGE2, GEDGE3, GEDGE4, GEDGE5, * GCORN0, GCORN1, GCORN2, GCORN3, GCORN4, GCORN5 C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C H2 = H*H F = H2*COEFU F2 = F*F A = -24.0E0 + 5.0E0*F - F2/4.0E0 B = 2.0E0 - F/24.0E0 + F2/48.0E0 C = 1.0E0 + 5.0E0*F/48.0E0 D = 0.0E0 BETA0 = 2.0E0 - F/4.0E0 BETA1 = 0.50E0 BETA2 = F/48.0E0 GAMMA0 = -12.0E0 + 11.0E0*F/12.0E0 GAMMA1 = F/12.0E0 GEDGE0 = -23.0E0/3.0E0 + 125.0E0*F/144.0E0 GEDGE1 = -6.0E0 + F/4.0E0 GEDGE2 = 2.0E0 - F/48.0E0 GEDGE3 = -1.0E0/3.0E0 - F/72.0E0 GEDGE4 = -0.50E0 + F/48.0E0 GEDGE5 = 0.50E0 + F/16.0E0 GCORN0 = (-2328.0E0 + 245.0E0*F)/144.0E0 GCORN1 = ( 141.0E0 - 14.0E0*F)/18.0E0 GCORN2 = ( -24.0E0 - F)/72.0E0 GCORN3 = ( 840.0E0 - 103.0E0*F)/36.0E0 GCORN4 = (-1752.0E0 + 221.0E0*F)/144.0E0 GCORN5 = ( 57.0E0 - 7.0E0*F)/9.0E0 C BETA0 = H2*BETA0 BETA1 = H2*BETA1 BETA2 = H2*BETA2 GAMMA0 = H*GAMMA0 GAMMA1 = H*GAMMA1 GEDGE0 = H*GEDGE0 GEDGE1 = H*GEDGE1 GEDGE2 = H*GEDGE2 GEDGE3 = H*GEDGE3 GEDGE4 = H*GEDGE4 GEDGE5 = H*GEDGE5 GCORN0 = H*GCORN0 GCORN1 = H*GCORN1 GCORN2 = H*GCORN2 GCORN3 = H*GCORN3 GCORN4 = H*GCORN4 GCORN5 = H*GCORN5 C LDXU = LMXU + 1 LDYU = LMYU + 1 NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 C HELMHZ = COEFU .NE. 0.0E0 HAVED = (BCTY(RIGHT) .EQ. DRCH) .OR. (BCTY(BOTTOM) .EQ. DRCH) .OR. * (BCTY(LEFT ) .EQ. DRCH) .OR. (BCTY(TOP ) .EQ. DRCH) .OR. * (BCTY(FRONT) .EQ. DRCH) .OR. (BCTY(BACK ) .EQ. DRCH) HAVEN = (BCTY(RIGHT) .EQ. NEUM) .OR. (BCTY(BOTTOM) .EQ. NEUM) .OR. * (BCTY(LEFT ) .EQ. NEUM) .OR. (BCTY(TOP ) .EQ. NEUM) .OR. * (BCTY(FRONT) .EQ. NEUM) .OR. (BCTY(BACK ) .EQ. NEUM) PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC PRDZ = BCTY(FRONT ) .EQ. PRDC C C C ----------------------------- C SHIFT VALUES OF G STORED IN U C ----------------------------- C DO 50 K=NZ,1,-1 DO 50 J=NY,1,-1 DO 50 I=NX,1,-1 U(I,J,K) = U(I-1,J-1,K-1) 50 CONTINUE C C C ----------------------------------------- C REFLECT FUNCTIONS G AND GH OUTSIDE DOMAIN C ----------------------------------------- C IF (HELMHZ) CALL REFL3(1,NX,NY,NZ,PRDX,PRDY,PRDZ,U,LMXU,LMYU) CALL REFL3(0,NXM1,NYM1,NZM1,PRDX,PRDY,PRDZ,GH,LMXGH,LMYGH) C C C ---------------------------- C DISCRETIZE RIGHT SIDE OF PDE C ---------------------------- C ISTRT = 2 ISTOP = NXM1 IF (BCTY(LEFT ) .NE. DRCH) ISTRT = 1 IF (BCTY(RIGHT ) .EQ. NEUM) ISTOP = NX JSTRT = 2 JSTOP = NYM1 IF (BCTY(BOTTOM) .NE. DRCH) JSTRT = 1 IF (BCTY(TOP ) .EQ. NEUM) JSTOP = NY KSTRT = 2 KSTOP = NZM1 IF (BCTY(BACK ) .NE. DRCH) KSTRT = 1 IF (BCTY(FRONT ) .EQ. NEUM) KSTOP = NZ C IF (HELMHZ) THEN C C CASE : HELMHOLTZ EQUATION C P0 = 0 DO 60 J=0,NY DO 60 I=0,NX WPLANE(I,J,P0) = U(I,J,KSTRT-1) 60 CONTINUE DO 100 K=KSTRT,KSTOP PM1 = P0 P0 = 1 - PM1 DO 70 J=0,NY DO 70 I=0,NX WPLANE(I,J,P0) = U(I,J,K) 70 CONTINUE DO 100 J=JSTRT,JSTOP DO 100 I=ISTRT,ISTOP C U(I,J,K) = BETA0*U(I,J,K) * + BETA1*( GH(I,J,K) + GH(I-1,J,K) + * GH(I,J-1,K) + GH(I-1,J,K-1) + * GH(I,J,K-1) + GH(I-1,J-1,K) + * GH(I,J-1,K-1) + GH(I-1,J-1,K-1) ) * + BETA2*( U(I+1,J,K) + WPLANE(I-1,J,P0) + * U(I,J+1,K) + WPLANE(I,J-1,P0) + * U(I,J,K+1) + WPLANE(I,J,PM1) ) 100 CONTINUE C ELSE C C CASE : POISSON EQUATION C DO 200 K=KSTRT,KSTOP DO 200 J=JSTRT,JSTOP DO 200 I=ISTRT,ISTOP C U(I,J,K) = BETA0*U(I,J,K) * + BETA1*( GH(I,J,K) + GH(I-1,J,K) + * GH(I,J-1,K) + GH(I-1,J,K-1) + * GH(I,J,K-1) + GH(I-1,J-1,K) + * GH(I,J-1,K-1) + GH(I-1,J-1,K-1) ) 200 CONTINUE C ENDIF C C C ------------------------- C REMOVE SHIFT FROM ARRAY U C ------------------------- C DO 250 K=0,NZM1 DO 250 J=0,NYM1 DO 250 I=0,NXM1 U(I,J,K) = U(I+1,J+1,K+1) 250 CONTINUE C C C ----------------------------------- C HANDLE POINTS ON NEUMANN BOUNDARIES C ----------------------------------- C IF (HAVEN) * CALL HD3N(NX,NY,NZ,BCTY,BD1,BD2,BD3,BD4,BD5,BD6,LDXBD,LDYBD, * U,LDXU,LDYU) C C C ---------------------------------------------- C ADJUST POINTS ADJACENT TO DIRICHLET BOUNDARIES C ---------------------------------------------- C IF (HAVED) THEN CALL STORD3(NX,NY,NZ,BCTY,BD1,BD2,BD3,BD4,BD5,BD6,LDXBD,LDYBD, * U,LDXU,LDYU) CALL FD3D(NX,NY,NZ,BCTY,B,C,D,U,LDXU,LDYU) ENDIF C C C ---- C EXIT C ---- C RETURN END SUBROUTINE HD3N (NX, NY, NZ, BCTY, BD1, BD2, BD3, BD4, BD5, BD6, * LDXBD, LDYBD, U, LDXU, LDYU) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C HD3N COMPUTES THE FOURTH ORDER FINITE DIFFERENCE DISCRETIZATION C AT ALL BOUNDARY POINTS OF A THREE-DIMENSIONAL RECTANGULAR DOMAIN C WHERE NEUMANN BOUNDARY CONDITIONS HAVE BEEN SPECIFIED. C C C P A R A M E T E R S C ------------------- C C NX,NY,NZ INTEGER SCALARS (INPUT) C SEE HFFT3A. C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C BD1, BD3 REAL ARRAYS OF SIZE LDYBD BY NZ (INPUT) C SEE HFFT3A. C C BD2, BD4 REAL ARRAYS OF SIZE LDXBD BY NZ (INPUT) C SEE HFFT3A. C C BD5, BD6 REAL ARRAYS OF SIZE LDXBD BY NY (INPUT) C SEE HFFT3A. C C LDXBD INTEGER SCALARS (INPUT) C LDYBD SEE HFFT3A. C C U REAL ARRAY OF SIZE LDXU BY LDYU BY NZ (OUTPUT) C ON EXIT, ENTRIES CORRESPONDING TO NEUMANN BOUNDARY C POINTS CONTAIN THE RIGHT HAND SIDE OF THE FINITE C DIFFERENCE DISCRETIZTION C C LDXU INTEGER SCALARS (INPUT) C LDYU SEE HFFT3A. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), NX, NY, NZ, LDXBD, LDYBD, LDXU, LDYU REAL * U(LDXU,LDYU,NZ), BD1(LDYBD,NZ), BD2(LDXBD,NZ), BD3(LDYBD,NZ), * BD4(LDXBD,NZ), BD5(LDXBD,NY), BD6(LDXBD,NY) C C ... LOCAL VARIABLES C LOGICAL PRDX, PRDY, PRDZ INTEGER I, J, K, NXM1, NYM1, NZM1, NXM2, NYM2, NZM2, * NXM3, NYM3, NZM3 REAL * GAMMA0, GAMMA1, * GEDGE0, GEDGE1, GEDGE2, GEDGE3, GEDGE4, GEDGE5, * GCORN0, GCORN1, GCORN2, GCORN3, GCORN4, GCORN5 C COMMON /HD3COM/ GAMMA0, GAMMA1, * GEDGE0, GEDGE1, GEDGE2, GEDGE3, GEDGE4, GEDGE5, * GCORN0, GCORN1, GCORN2, GCORN3, GCORN4, GCORN5 C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C --------------- C INITIALIZATIONS C --------------- C NXM1 = NX - 1 NYM1 = NY - 1 NZM1 = NZ - 1 NXM2 = NX - 2 NYM2 = NY - 2 NZM2 = NZ - 2 NXM3 = NX - 3 NYM3 = NY - 3 NZM3 = NZ - 3 PRDX = BCTY(RIGHT ) .EQ. PRDC PRDY = BCTY(BOTTOM) .EQ. PRDC PRDZ = BCTY(FRONT ) .EQ. PRDC C C C ------------------------------ C HANDLE POINTS ON NEUMANN SIDES C ------------------------------ C C ... RIGHT SIDE C IF (BCTY(RIGHT) .EQ. NEUM) THEN DO 205 K=2,NZM1 DO 205 J=2,NYM1 U(NX,J,K) = U(NX,J,K) * + GAMMA0*BD1(J,K) * + GAMMA1*(BD1(J+1,K) + BD1(J-1,K) + * BD1(J,K+1) + BD1(J,K-1) ) 205 CONTINUE ENDIF C C ... LEFT SIDE C IF (BCTY(LEFT) .EQ. NEUM) THEN DO 215 K=2,NZM1 DO 215 J=2,NYM1 U(1,J,K) = U(1,J,K) * - GAMMA0*BD3(J,K) * - GAMMA1*(BD3(J+1,K) + BD3(J-1,K) + * BD3(J,K+1) + BD3(J,K-1) ) 215 CONTINUE ENDIF C C ... TOP SIDE C IF (BCTY(TOP) .EQ. NEUM) THEN DO 225 K=2,NZM1 DO 225 I=2,NXM1 U(I,NY,K) = U(I,NY,K) * + GAMMA0*BD4(I,K) * + GAMMA1*(BD4(I+1,K) + BD4(I-1,K) + * BD4(I,K+1) + BD4(I,K-1) ) 225 CONTINUE ENDIF C C ... BOTTOM SIDE C IF (BCTY(BOTTOM) .EQ. NEUM) THEN DO 235 K=2,NZM1 DO 235 I=2,NXM1 U(I,1,K) = U(I,1,K) * - GAMMA0*BD2(I,K) * - GAMMA1*(BD2(I+1,K) + BD2(I-1,K) + * BD2(I,K+1) + BD2(I,K-1) ) 235 CONTINUE ENDIF C C ... FRONT SIDE C IF (BCTY(FRONT) .EQ. NEUM) THEN DO 245 J=2,NYM1 DO 245 I=2,NXM1 U(I,J,NZ) = U(I,J,NZ) * + GAMMA0*BD5(I,J) * + GAMMA1*(BD5(I+1,J) + BD5(I-1,J) + * BD5(I,J+1) + BD5(I,J-1) ) 245 CONTINUE ENDIF C C ... BACK SIDE C IF (BCTY(BACK) .EQ. NEUM) THEN DO 255 J=2,NYM1 DO 255 I=2,NXM1 U(I,J,1) = U(I,J,1) * - GAMMA0*BD6(I,J) * - GAMMA1*(BD6(I+1,J) + BD6(I-1,J) + * BD6(I,J+1) + BD6(I,J-1) ) 255 CONTINUE ENDIF C C C ------------------------------ C HANDLE POINTS ON NEUMANN EDGES C ------------------------------ C C ... RIGHT TOP EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 265 K=2,NZM1 U(NX,NY,K) = U(NX,NY,K) * + GEDGE0*( BD1(NY,K) + BD4(NX,K) ) * + GEDGE1*( BD1(NYM1,K) + BD4(NXM1,K) ) * + GEDGE2*( BD1(NYM2,K) + BD4(NXM2,K) ) * + GEDGE3*( BD1(NYM3,K) + BD4(NXM3,K) ) * + GEDGE4*( BD1(NY,K+1) + BD4(NX,K+1) * + BD1(NY,K-1) + BD4(NX,K-1) ) * + GEDGE5*( BD1(NYM1,K+1) + BD4(NXM1,K+1) * + BD1(NYM1,K-1) + BD4(NXM1,K-1) ) 265 CONTINUE ENDIF C C ... RIGHT BOTTOM EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) THEN DO 275 K=2,NZM1 U(NX,1,K) = U(NX,1,K) * + GEDGE0*( BD1(1,K) - BD2(NX,K) ) * + GEDGE1*( BD1(2,K) - BD2(NXM1,K) ) * + GEDGE2*( BD1(3,K) - BD2(NXM2,K) ) * + GEDGE3*( BD1(4,K) - BD2(NXM3,K) ) * + GEDGE4*( BD1(1,K+1) - BD2(NX,K+1) * + BD1(1,K-1) - BD2(NX,K-1) ) * + GEDGE5*( BD1(2,K+1) - BD2(NXM1,K+1) * + BD1(2,K-1) - BD2(NXM1,K-1) ) 275 CONTINUE ENDIF C C ... LEFT TOP EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 285 K=2,NZM1 U(1,NY,K) = U(1,NY,K) * + GEDGE0*( -BD3(NY,K) + BD4(1,K) ) * + GEDGE1*( -BD3(NYM1,K) + BD4(2,K) ) * + GEDGE2*( -BD3(NYM2,K) + BD4(3,K) ) * + GEDGE3*( -BD3(NYM3,K) + BD4(4,K) ) * + GEDGE4*( -BD3(NY,K+1) + BD4(1,K+1) * -BD3(NY,K-1) + BD4(1,K-1) ) * + GEDGE5*( -BD3(NYM1,K+1) + BD4(2,K+1) * -BD3(NYM1,K-1) + BD4(2,K-1) ) 285 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(BOTTOM) .EQ. NEUM)) THEN DO 295 K=2,NZM1 U(1,1,K) = U(1,1,K) * + GEDGE0*( -BD3(1,K) - BD2(1,K) ) * + GEDGE1*( -BD3(2,K) - BD2(2,K) ) * + GEDGE2*( -BD3(3,K) - BD2(3,K) ) * + GEDGE3*( -BD3(4,K) - BD2(4,K) ) * + GEDGE4*( -BD3(1,K+1) - BD2(1,K+1) * -BD3(1,K-1) - BD2(1,K-1) ) * + GEDGE5*( -BD3(2,K+1) - BD2(2,K+1) * -BD3(2,K-1) - BD2(2,K-1) ) 295 CONTINUE ENDIF C C ... RIGHT FRONT EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 305 J=2,NYM1 U(NX,J,NZ) = U(NX,J,NZ) * + GEDGE0*( BD1(J,NZ) + BD5(NX,J) ) * + GEDGE1*( BD1(J,NZM1) + BD5(NXM1,J) ) * + GEDGE2*( BD1(J,NZM2) + BD5(NXM2,J) ) * + GEDGE3*( BD1(J,NZM3) + BD5(NXM3,J) ) * + GEDGE4*( BD1(J+1,NZ) + BD5(NX,J+1) * + BD1(J-1,NZ) + BD5(NX,J-1) ) * + GEDGE5*( BD1(J+1,NZM1) + BD5(NXM1,J+1) * + BD1(J-1,NZM1) + BD5(NXM1,J-1) ) 305 CONTINUE ENDIF C C ... RIGHT BACK EDGE C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 315 J=2,NYM1 U(NX,J,1) = U(NX,J,1) * + GEDGE0*( BD1(J,1) - BD6(NX,J) ) * + GEDGE1*( BD1(J,2) - BD6(NXM1,J) ) * + GEDGE2*( BD1(J,3) - BD6(NXM2,J) ) * + GEDGE3*( BD1(J,4) - BD6(NXM3,J) ) * + GEDGE4*( BD1(J+1,1) - BD6(NX,J+1) * + BD1(J-1,1) - BD6(NX,J-1) ) * + GEDGE5*( BD1(J+1,2) - BD6(NXM1,J+1) * + BD1(J-1,2) - BD6(NXM1,J-1) ) 315 CONTINUE ENDIF C C ... LEFT FRONT EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 325 J=2,NYM1 U(1,J,NZ) = U(1,J,NZ) * + GEDGE0*( -BD3(J,NZ) + BD5(1,J) ) * + GEDGE1*( -BD3(J,NZM1) + BD5(2,J) ) * + GEDGE2*( -BD3(J,NZM2) + BD5(3,J) ) * + GEDGE3*( -BD3(J,NZM3) + BD5(4,J) ) * + GEDGE4*( -BD3(J+1,NZ) + BD5(1,J+1) * -BD3(J-1,NZ) + BD5(1,J-1) ) * + GEDGE5*( -BD3(J+1,NZM1) + BD5(2,J+1) * -BD3(J-1,NZM1) + BD5(2,J-1) ) 325 CONTINUE ENDIF C C ... LEFT FRON EDGE C IF ((BCTY(LEFT) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 335 J=2,NYM1 U(1,J,1) = U(1,J,1) * + GEDGE0*( -BD3(J,1) - BD6(1,J) ) * + GEDGE1*( -BD3(J,2) - BD6(2,J) ) * + GEDGE2*( -BD3(J,3) - BD6(3,J) ) * + GEDGE3*( -BD3(J,4) - BD6(4,J) ) * + GEDGE4*( -BD3(J+1,1) - BD6(1,J+1) * -BD3(J-1,1) - BD6(1,J-1) ) * + GEDGE5*( -BD3(J+1,2) - BD6(2,J+1) * -BD3(J-1,2) - BD6(2,J-1) ) 335 CONTINUE ENDIF C C ... TOP FRONT EDGE C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 345 I=2,NXM1 U(I,NY,NZ) = U(I,NY,NZ) * + GEDGE0*( BD4(I,NZ) + BD5(I,NY) ) * + GEDGE1*( BD4(I,NZM1) + BD5(I,NYM1) ) * + GEDGE2*( BD4(I,NZM2) + BD5(I,NYM2) ) * + GEDGE3*( BD4(I,NZM3) + BD5(I,NYM3) ) * + GEDGE4*( BD4(I+1,NZ) + BD5(I+1,NY) * + BD4(I-1,NZ) + BD5(I-1,NY) ) * + GEDGE5*( BD4(I+1,NZM1) + BD5(I+1,NYM1) * + BD4(I-1,NZM1) + BD5(I-1,NYM1) ) 345 CONTINUE ENDIF C C ... TOP BACK EDGE C IF ((BCTY(TOP) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 355 I=2,NXM1 U(I,NY,1) = U(I,NY,1) * + GEDGE0*( BD4(I,1) - BD6(I,NY) ) * + GEDGE1*( BD4(I,2) - BD6(I,NYM1) ) * + GEDGE2*( BD4(I,3) - BD6(I,NYM2) ) * + GEDGE3*( BD4(I,4) - BD6(I,NYM3) ) * + GEDGE4*( BD4(I+1,1) - BD6(I+1,NY) * + BD4(I-1,1) - BD6(I-1,NY) ) * + GEDGE5*( BD4(I+1,2) - BD6(I+1,NYM1) * + BD4(I-1,2) - BD6(I-1,NYM1) ) 355 CONTINUE ENDIF C C ... BOTTOM FRONT EDGE C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 365 I=2,NXM1 U(I,1,NZ) = U(I,1,NZ) * + GEDGE0*( -BD2(I,NZ) + BD5(I,1) ) * + GEDGE1*( -BD2(I,NZM1) + BD5(I,2) ) * + GEDGE2*( -BD2(I,NZM2) + BD5(I,3) ) * + GEDGE3*( -BD2(I,NZM3) + BD5(I,4) ) * + GEDGE4*( -BD2(I+1,NZ) + BD5(I+1,1) * -BD2(I-1,NZ) + BD5(I-1,1) ) * + GEDGE5*( -BD2(I+1,NZM1) + BD5(I+1,2) * -BD2(I-1,NZM1) + BD5(I-1,2) ) 365 CONTINUE ENDIF C C ... BOTTOM BACK EDGE C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 375 I=2,NXM1 U(I,1,1) = U(I,1,1) * + GEDGE0*( -BD2(I,1) - BD6(I,1) ) * + GEDGE1*( -BD2(I,2) - BD6(I,2) ) * + GEDGE2*( -BD2(I,3) - BD6(I,3) ) * + GEDGE3*( -BD2(I,4) - BD6(I,4) ) * + GEDGE4*( -BD2(I+1,1) - BD6(I+1,1) * -BD2(I-1,1) - BD6(I-1,1) ) * + GEDGE5*( -BD2(I+1,2) - BD6(I+1,2) * -BD2(I-1,2) - BD6(I-1,2) ) 375 CONTINUE ENDIF C C C -------------------------------- C HANDLE POINTS AT NEUMANN CORNERS C -------------------------------- C C ... RIGHT TOP FRONT CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(NX,NY,NZ) = U(NX,NY,NZ) * + GCORN0*( BD1(NYM1,NZ) + BD1(NY,NZM1) * + BD4(NXM1,NZ) + BD4(NX,NZM1) * + BD5(NXM1,NY) + BD5(NX,NYM1) ) * + GCORN1*( BD1(NYM2,NZ) + BD1(NY,NZM2) * + BD4(NXM2,NZ) + BD4(NX,NZM2) * + BD5(NXM2,NY) + BD5(NX,NYM2) ) * + GCORN2*( BD1(NYM3,NZ) + BD1(NY,NZM3) * + BD4(NXM3,NZ) + BD4(NX,NZM3) * + BD5(NXM3,NY) + BD5(NX,NYM3) ) * + GCORN3*( BD1(NYM1,NZM1)+BD4(NXM1,NZM1)+BD5(NXM1,NYM1)) * + GCORN4*( BD1(NYM1,NZM2) + BD1(NYM2,NZM1) * + BD4(NXM1,NZM2) + BD4(NXM2,NZM1) * + BD5(NXM1,NYM2) + BD5(NXM2,NYM1) ) * + GCORN5*( BD1(NYM2,NZM2)+BD4(NXM2,NZM2)+BD5(NXM2,NYM2)) C C ... LEFT TOP FRONT CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM)) * U(1,NY,NZ) = U(1,NY,NZ) * + GCORN0*(-BD3(NYM1,NZ) - BD3(NY,NZM1) * + BD4(2,NZ) + BD4(1,NZM1) * + BD5(2,NY) + BD5(1,NYM1) ) * + GCORN1*(-BD3(NYM2,NZ) - BD3(NY,NZM2) * + BD4(3,NZ) + BD4(1,NZM2) * + BD5(3,NY) + BD5(1,NYM2) ) * + GCORN2*(-BD3(NYM3,NZ) - BD3(NY,NZM3) * + BD4(4,NZ) + BD4(1,NZM3) * + BD5(4,NY) + BD5(1,NYM3) ) * + GCORN3*(-BD3(NYM1,NZM1)+BD4(2,NZM1)+BD5(2,NYM1)) * + GCORN4*(-BD3(NYM1,NZM2) - BD3(NYM2,NZM1) * + BD4(2,NZM2) + BD4(3,NZM1) * + BD5(2,NYM2) + BD5(3,NYM1) ) * + GCORN5*(-BD3(NYM2,NZM2)+BD4(3,NZM2)+BD5(3,NYM2)) C C ... RIGHT BOTTOM FRONT CORNER C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM)) * U(NX,1,NZ) = U(NX,1,NZ) * + GCORN0*( BD1(2,NZ) + BD1(1,NZM1) * - BD2(NXM1,NZ) - BD2(NX,NZM1) * + BD5(NXM1,1) + BD5(NX,2) ) * + GCORN1*( BD1(3,NZ) + BD1(1,NZM2) * - BD2(NXM2,NZ) - BD2(NX,NZM2) * + BD5(NXM2,1) + BD5(NX,3) ) * + GCORN2*( BD1(4,NZ) + BD1(1,NZM3) * - BD2(NXM3,NZ) - BD2(NX,NZM3) * + BD5(NXM3,1) + BD5(NX,4) ) * + GCORN3*( BD1(2,NZM1)-BD2(NXM1,NZM1)+BD5(NXM1,2)) * + GCORN4*( BD1(2,NZM2) + BD1(3,NZM1) * - BD2(NXM1,NZM2) - BD2(NXM2,NZM1) * + BD5(NXM1,3) + BD5(NXM2,2) ) * + GCORN5*( BD1(3,NZM2)-BD2(NXM2,NZM2)+BD5(NXM2,3)) C C ... LEFT BOTTOM FRONT CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM)) * U(1,1,NZ) = U(1,1,NZ) * + GCORN0*(-BD3(2,NZ) - BD3(1,NZM1) * - BD2(2,NZ) - BD2(1,NZM1) * + BD5(2,1) + BD5(1,2) ) * + GCORN1*(-BD3(3,NZ) - BD3(1,NZM2) * - BD2(3,NZ) - BD2(1,NZM2) * + BD5(3,1) + BD5(1,3) ) * + GCORN2*(-BD3(4,NZ) - BD3(1,NZM3) * - BD2(4,NZ) - BD2(1,NZM3) * + BD5(4,1) + BD5(1,4) ) * + GCORN3*(-BD3(2,NZM1)-BD2(2,NZM1)+BD5(2,2)) * + GCORN4*(-BD3(2,NZM2) - BD3(3,NZM1) * - BD2(2,NZM2) - BD2(3,NZM1) * + BD5(2,3) + BD5(3,2) ) * + GCORN5*(-BD3(3,NZM2)-BD2(3,NZM2)+BD5(3,3)) C C ... RIGHT TOP BACK CORNER C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(NX,NY,1) = U(NX,NY,1) * + GCORN0*( BD1(NYM1,1) + BD1(NY,2) * + BD4(NXM1,1) + BD4(NX,2) * - BD6(NXM1,NY) - BD6(NX,NYM1) ) * + GCORN1*( BD1(NYM2,1) + BD1(NY,3) * + BD4(NXM2,1) + BD4(NX,3) * - BD6(NXM2,NY) - BD6(NX,NYM2) ) * + GCORN2*( BD1(NYM3,1) + BD1(NY,4) * + BD4(NXM3,1) + BD4(NX,4) * - BD6(NXM3,NY) - BD6(NX,NYM3) ) * + GCORN3*( BD1(NYM1,2)+BD4(NXM1,2)-BD6(NXM1,NYM1)) * + GCORN4*( BD1(NYM1,3) + BD1(NYM2,2) * + BD4(NXM1,3) + BD4(NXM2,2) * - BD6(NXM1,NYM2) - BD6(NXM2,NYM1) ) * + GCORN5*( BD1(NYM2,3)+BD4(NXM2,3)-BD6(NXM2,NYM2)) C C ... LEFT TOP BACK CORNER C IF ((BCTY(LEFT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(BACK) .EQ. NEUM)) * U(1,NY,1) = U(1,NY,1) * + GCORN0*(-BD3(NYM1,1) - BD3(NY,2) * + BD4(2,1) + BD4(1,2) * - BD6(2,NY) - BD6(1,NYM1) ) * + GCORN1*(-BD3(NYM2,1) - BD3(NY,3) * + BD4(3,1) + BD4(1,3) * - BD6(3,NY) - BD6(1,NYM2) ) * + GCORN2*(-BD3(NYM3,1) - BD3(NY,4) * + BD4(4,1) + BD4(1,4) * - BD6(4,NY) - BD6(1,NYM3) ) * + GCORN3*(-BD3(NYM1,2)+BD4(2,2)-BD6(2,NYM1)) * + GCORN4*(-BD3(NYM1,3) - BD3(NYM2,2) * + BD4(2,3) + BD4(3,2) * - BD6(2,NYM2) - BD6(3,NYM1) ) * + GCORN5*(-BD3(NYM2,3)+BD4(3,3)-BD6(3,NYM2)) C C ... RIGHT BOTTOM BACK CORNER C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(NX,1,1) = U(NX,1,1) * + GCORN0*( BD1(2,1) + BD1(1,2) * - BD2(NXM1,1) - BD2(NX,2) * - BD6(NXM1,1) - BD6(NX,2) ) * + GCORN1*( BD1(3,1) + BD1(1,3) * - BD2(NXM2,1) - BD2(NX,3) * - BD6(NXM2,1) - BD6(NX,3) ) * + GCORN2*( BD1(4,1) + BD1(1,4) * - BD2(NXM3,1) - BD2(NX,4) * - BD6(NXM3,1) - BD6(NX,4) ) * + GCORN3*( BD1(2,2)-BD2(NXM1,2)-BD6(NXM1,2)) * + GCORN4*( BD1(2,3) + BD1(3,2) * - BD2(NXM1,3) - BD2(NXM2,2) * - BD6(NXM1,3) - BD6(NXM2,2) ) * + GCORN5*( BD1(3,3)-BD2(NXM2,3)-BD6(NXM2,3)) C C ... LEFT BOTTOM BACK CORNER C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK ) .EQ. NEUM)) * U(1,1,1) = U(1,1,1) * + GCORN0*(-BD3(2,1) - BD3(1,2) * - BD2(2,1) - BD2(1,2) * - BD6(2,1) - BD6(1,2) ) * + GCORN1*(-BD3(3,1) - BD3(1,3) * - BD2(3,1) - BD2(1,3) * - BD6(3,1) - BD6(1,3) ) * + GCORN2*(-BD3(4,1) - BD3(1,4) * - BD2(4,1) - BD2(1,4) * - BD6(4,1) - BD6(1,4) ) * + GCORN3*(-BD3(2,2)-BD2(2,2)-BD6(2,2)) * + GCORN4*(-BD3(2,3) - BD3(3,2) * - BD2(2,3) - BD2(3,2) * - BD6(2,3) - BD6(3,2) ) * + GCORN5*(-BD3(3,3)-BD2(3,3)-BD6(3,3)) C C IF (.NOT. (PRDX .OR. PRDY .OR. PRDZ)) GO TO 999 C C C --------------------------------------- C HANDLE POINTS ON NEUMANN/PERIODIC EDGES C --------------------------------------- C C ... BOTTOM BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) THEN DO 405 I=2,NXM1 U(I,1,1) = U(I,1,1) * - GAMMA0*BD2(I,1) * - GAMMA1*(BD2(I+1,1) + BD2(I-1,1) + * BD2(I,2) + BD2(I,NZM1) ) 405 CONTINUE ENDIF C C ... BOTTOM BACK EDGE (PERIODIC/NEUMANN) C IF (PRDY .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 415 I=2,NXM1 U(I,1,1) = U(I,1,1) * - GAMMA0*BD6(I,1) * - GAMMA1*(BD6(I+1,1) + BD6(I-1,1) + * BD6(I,2) + BD6(I,NYM1) ) 415 CONTINUE ENDIF C C ... LEFT BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDZ) THEN DO 425 J=2,NYM1 U(1,J,1) = U(1,J,1) * - GAMMA0*BD3(J,1) * - GAMMA1*(BD3(J+1,1) + BD3(J-1,1) + * BD3(J,2) + BD3(J,NZM1) ) 425 CONTINUE ENDIF C C ... LEFT BACK EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(BACK) .EQ. NEUM)) THEN DO 435 J=2,NYM1 U(1,J,1) = U(1,J,1) * - GAMMA0*BD6(1,J) * - GAMMA1*(BD6(2,J) + BD6(NXM1,J) + * BD6(1,J+1) + BD6(1,J-1) ) 435 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE (NEUMANN/PERIODIC) C IF ((BCTY(BOTTOM) .EQ. NEUM) .AND. PRDX) THEN DO 445 K=2,NZM1 U(1,1,K) = U(1,1,K) * - GAMMA0*BD2(1,K) * - GAMMA1*(BD2(2,K) + BD2(NXM1,K) + * BD2(1,K+1) + BD2(1,K-1) ) 445 CONTINUE ENDIF C C ... LEFT BOTTOM EDGE (NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY) THEN DO 455 K=2,NZM1 U(1,1,K) = U(1,1,K) * - GAMMA0*BD3(1,K) * - GAMMA1*(BD3(2,K) + BD3(NYM1,K) + * BD3(1,K+1) + BD3(1,K-1) ) 455 CONTINUE ENDIF C C ... LEFT TOP EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM)) THEN DO 465 K=2,NZM1 U(1,NY,K) = U(1,NY,K) * + GAMMA0*BD4(1,K) * + GAMMA1*(BD4(2,K) + BD4(NXM1,K) + * BD4(1,K+1) + BD4(1,K-1) ) 465 CONTINUE ENDIF C C ... LEFT FRONT EDGE (PERIODIC/NEUMANN) C IF (PRDX .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 475 J=2,NYM1 U(1,J,NZ) = U(1,J,NZ) * + GAMMA0*BD5(1,J) * + GAMMA1*(BD5(2,J) + BD5(NXM1,J) + * BD5(1,J+1) + BD5(1,J-1) ) 475 CONTINUE ENDIF C C ... RIGHT BOTTOM EDGE (NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY) THEN DO 485 K=2,NZM1 U(NX,1,K) = U(NX,1,K) * + GAMMA0*BD1(1,K) * + GAMMA1*(BD1(2,K) + BD1(NYM1,K) + * BD1(1,K+1) + BD1(1,K-1) ) 485 CONTINUE ENDIF C C ... BOTTOM FRONT EDGE (PERIODIC/NEUMANN) C IF (PRDY .AND. (BCTY(FRONT) .EQ. NEUM)) THEN DO 495 I=2,NXM1 U(I,1,NZ) = U(I,1,NZ) * + GAMMA0*BD5(I,1) * + GAMMA1*(BD5(I+1,1) + BD5(I-1,1) + * BD5(I,2) + BD5(I,NYM1) ) 495 CONTINUE ENDIF C C ... RIGHT BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDZ) THEN DO 505 J=2,NYM1 U(NX,J,1) = U(NX,J,1) * + GAMMA0*BD1(J,1) * + GAMMA1*(BD1(J+1,1) + BD1(J-1,1) + * BD1(J,2) + BD1(J,NZM1) ) 505 CONTINUE ENDIF C C ... TOP BACK EDGE (NEUMANN/PERIODIC) C IF ((BCTY(TOP) .EQ. NEUM) .AND. PRDZ) THEN DO 515 I=2,NXM1 U(I,NY,1) = U(I,NY,1) * + GAMMA0*BD4(I,1) * + GAMMA1*(BD4(I+1,1) + BD4(I-1,1) + * BD4(I,2) + BD4(I,NZM1) ) 515 CONTINUE ENDIF C C C ----------------------------------------- C HANDLE NEUMANN/POINTS AT PERIODIC CORNERS C ----------------------------------------- C C ... LEFT BOTTOM BACK CORNER (NEUMANN/PERIODIC/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY .AND. PRDZ) THEN U(1,1,1) = U(1,1,1) * - GAMMA0*BD3(1,1) * - GAMMA1*(BD3(2,1) + BD3(NYM1,1) + * BD3(1,2) + BD3(1,NZM1) ) ENDIF C C ... LEFT BOTTOM BACK CORNER (PERIODIC/NEUMANN/PERIODIC) C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) THEN U(1,1,1) = U(1,1,1) * - GAMMA0*BD2(1,1) * - GAMMA1*(BD2(2,1) + BD2(NXM1,1) + * BD2(1,2) + BD2(1,NZM1) ) ENDIF C C ... LEFT BOTTOM BACK CORNER (PERIODIC/PERIODIC/NEUMANN) C IF (PRDX .AND. PRDY .AND. (BCTY(BACK) .EQ. NEUM)) THEN U(1,1,1) = U(1,1,1) * - GAMMA0*BD6(1,1) * - GAMMA1*(BD6(2,1) + BD6(NXM1,1) + * BD6(1,2) + BD6(1,NYM1) ) ENDIF C C ... LEFT BOTTOM BACK CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(BACK) .EQ. NEUM) ) THEN U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD2(1,1) - BD6(1,1) ) * + GEDGE1*( -BD2(1,2) - BD6(1,2) ) * + GEDGE2*( -BD2(1,3) - BD6(1,3) ) * + GEDGE3*( -BD2(1,4) - BD6(1,4) ) * + GEDGE4*( -BD2(2,1) - BD6(2,1) * -BD2(NXM1,1) - BD6(NXM1,1) ) * + GEDGE5*( -BD2(2,2) - BD6(2,2) * -BD2(NXM1,2) - BD6(NXM1,2) ) ENDIF C C ... LEFT BOTTOM BACK CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(BACK) .EQ. NEUM) ) THEN U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD3(1,1) - BD6(1,1) ) * + GEDGE1*( -BD3(1,2) - BD6(2,1) ) * + GEDGE2*( -BD3(1,3) - BD6(3,1) ) * + GEDGE3*( -BD3(1,4) - BD6(4,1) ) * + GEDGE4*( -BD3(2,1) - BD6(1,2) * -BD3(NYM1,1) - BD6(1,NYM1) ) * + GEDGE5*( -BD3(2,2) - BD6(2,2) * -BD3(NYM1,2) - BD6(2,NYM1) ) ENDIF C C ... LEFT BOTTOM BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) THEN U(1,1,1) = U(1,1,1) * + GEDGE0*( -BD3(1,1) - BD2(1,1) ) * + GEDGE1*( -BD3(2,1) - BD2(2,1) ) * + GEDGE2*( -BD3(3,1) - BD2(3,1) ) * + GEDGE3*( -BD3(4,1) - BD2(4,1) ) * + GEDGE4*( -BD3(1,2) - BD2(1,2) * -BD3(1,NZM1) - BD2(1,NZM1) ) * + GEDGE5*( -BD3(2,2) - BD2(2,2) * -BD3(2,NZM1) - BD2(2,NZM1) ) ENDIF C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/PERIODIC/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. PRDZ) THEN U(NX,1,1) = U(NX,1,1) * + GAMMA0*BD1(1,1) * + GAMMA1*(BD1(2,1) + BD1(NYM1,1) + * BD1(1,2) + BD1(1,NZM1) ) ENDIF C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(RIGHT ) .EQ. NEUM) .AND. * (BCTY(BOTTOM) .EQ. NEUM) .AND. PRDZ) THEN U(NX,1,1) = U(NX,1,1) * + GEDGE0*( BD1(1,1) - BD2(NX,1) ) * + GEDGE1*( BD1(2,1) - BD2(NXM1,1) ) * + GEDGE2*( BD1(3,1) - BD2(NXM2,1) ) * + GEDGE3*( BD1(4,1) - BD2(NXM3,1) ) * + GEDGE4*( BD1(1,2) - BD2(NX,2) * + BD1(1,NZM1) - BD2(NX,NZM1) ) * + GEDGE5*( BD1(2,2) - BD2(NXM1,2) * + BD1(2,NZM1) - BD2(NXM1,NZM1) ) ENDIF C C ... RIGHT BOTTOM BACK CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(BACK ) .EQ. NEUM) ) THEN U(NX,1,1) = U(NX,1,1) * + GEDGE0*( BD1(1,1) - BD6(NX,1) ) * + GEDGE1*( BD1(1,2) - BD6(NXM1,1) ) * + GEDGE2*( BD1(1,3) - BD6(NXM2,1) ) * + GEDGE3*( BD1(1,4) - BD6(NXM3,1) ) * + GEDGE4*( BD1(2,1) - BD6(NX,2) * + BD1(NYM1,1) - BD6(NX,NYM1) ) * + GEDGE5*( BD1(2,2) - BD6(NXM1,2) * + BD1(NYM1,2) - BD6(NXM1,NYM1) ) ENDIF C C ... LEFT TOP BACK CORNER (PERIODIC/NEUMANN/PERIODIC) C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM) .AND. PRDZ) THEN U(1,NY,1) = U(1,NY,1) * + GAMMA0*BD4(1,1) * + GAMMA1*(BD4(2,1) + BD4(NXM1,1) + * BD4(1,2) + BD4(1,NZM1) ) ENDIF C C ... LEFT TOP BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(LEFT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. PRDZ) THEN U(1,NY,1) = U(1,NY,1) * + GEDGE0*( -BD3(NY,1) + BD4(1,1) ) * + GEDGE1*( -BD3(NYM1,1) + BD4(2,1) ) * + GEDGE2*( -BD3(NYM2,1) + BD4(3,1) ) * + GEDGE3*( -BD3(NYM3,1) + BD4(4,1) ) * + GEDGE4*( -BD3(NY,2) + BD4(1,2) * -BD3(NY,NZM1) + BD4(1,NZM1) ) * + GEDGE5*( -BD3(NYM1,2) + BD4(2,2) * -BD3(NYM1,NZM1) + BD4(2,NZM1) ) ENDIF C C ... LEFT TOP BACK CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. (BCTY(TOP) .EQ. NEUM) .AND. * (BCTY(BACK) .EQ. NEUM) ) THEN U(1,NY,1) = U(1,NY,1) * + GEDGE0*( BD4(1,1) - BD6(1,NY) ) * + GEDGE1*( BD4(1,2) - BD6(1,NYM1) ) * + GEDGE2*( BD4(1,3) - BD6(1,NYM2) ) * + GEDGE3*( BD4(1,4) - BD6(1,NYM3) ) * + GEDGE4*( BD4(2,1) - BD6(2,NY) * + BD4(NXM1,1) - BD6(NXM1,NY) ) * + GEDGE5*( BD4(2,2) - BD6(2,NYM1) * + BD4(NXM1,2) - BD6(NXM1,NYM1) ) ENDIF C C ... LEFT BOTTOM FRONT CORNER (PERIODIC/PERIODIC/NEUMANN) C IF (PRDX .AND. PRDY .AND. (BCTY(FRONT) .EQ. NEUM)) THEN U(1,1,NZ) = U(1,1,NZ) * + GAMMA0*BD5(1,1) * + GAMMA1*(BD5(2,1) + BD5(NXM1,1) + * BD5(1,2) + BD5(1,NYM1) ) ENDIF C C ... LEFT BOTTOM FRONT CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(LEFT ) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(FRONT) .EQ. NEUM) ) THEN U(1,1,NZ) = U(1,1,NZ) * + GEDGE0*( -BD3(1,NZ) + BD5(1,1) ) * + GEDGE1*( -BD3(1,NZM1) + BD5(2,1) ) * + GEDGE2*( -BD3(1,NZM2) + BD5(3,1) ) * + GEDGE3*( -BD3(1,NZM3) + BD5(4,1) ) * + GEDGE4*( -BD3(2,NZ) + BD5(1,2) * -BD3(NYM1,NZ) + BD5(1,NYM1) ) * + GEDGE5*( -BD3(2,NZM1) + BD5(2,2) * -BD3(NYM1,NZM1) + BD5(2,NYM1) ) ENDIF C C ... LEFT BOTTOM FRONT CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. (BCTY(BOTTOM) .EQ. NEUM) .AND. * (BCTY(FRONT ) .EQ. NEUM) ) THEN U(1,1,NZ) = U(1,1,NZ) * + GEDGE0*( -BD2(1,NZ) + BD5(1,1) ) * + GEDGE1*( -BD2(1,NZM1) + BD5(1,2) ) * + GEDGE2*( -BD2(1,NZM2) + BD5(1,3) ) * + GEDGE3*( -BD2(1,NZM3) + BD5(1,4) ) * + GEDGE4*( -BD2(2,NZ) + BD5(2,1) * -BD2(NXM1,NZ) + BD5(NXM1,1) ) * + GEDGE5*( -BD2(2,NZM1) + BD5(2,2) * -BD2(NXM1,NZM1) + BD5(NXM1,2) ) ENDIF C C ... LEFT TOP FRONT CORNER (PERIODIC/NEUMANN/NEUMANN) C IF (PRDX .AND. (BCTY(TOP ) .EQ. NEUM) .AND. * (BCTY(FRONT) .EQ. NEUM) ) THEN U(1,NY,NZ) = U(1,NY,NZ) * + GEDGE0*( BD4(1,NZ) + BD5(1,NY) ) * + GEDGE1*( BD4(1,NZM1) + BD5(1,NYM1) ) * + GEDGE2*( BD4(1,NZM2) + BD5(1,NYM2) ) * + GEDGE3*( BD4(1,NZM3) + BD5(1,NYM3) ) * + GEDGE4*( BD4(2,NZ) + BD5(2,NY) * + BD4(NXM1,NZ) + BD5(NXM1,NY) ) * + GEDGE5*( BD4(2,NZM1) + BD5(2,NYM1) * + BD4(NXM1,NZM1) + BD5(NXM1,NYM1) ) ENDIF C C ... RIGHT BOTTOM FRONT CORNER (NEUMANN/PERIODIC/NEUMANN) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. PRDY .AND. * (BCTY(FRONT) .EQ. NEUM) ) THEN U(NX,1,NZ) = U(NX,1,NZ) * + GEDGE0*( BD1(1,NZ) + BD5(NX,1) ) * + GEDGE1*( BD1(1,NZM1) + BD5(NXM1,1) ) * + GEDGE2*( BD1(1,NZM2) + BD5(NXM2,1) ) * + GEDGE3*( BD1(1,NZM3) + BD5(NXM3,1) ) * + GEDGE4*( BD1(2,NZ) + BD5(NX,2) * + BD1(NYM1,NZ) + BD5(NX,NYM1) ) * + GEDGE5*( BD1(2,NZM1) + BD5(NXM1,2) * + BD1(NYM1,NZM1) + BD5(NXM1,NYM1) ) ENDIF C C ... RIGHT TOP BACK CORNER (NEUMANN/NEUMANN/PERIODIC) C IF ((BCTY(RIGHT) .EQ. NEUM) .AND. * (BCTY(TOP ) .EQ. NEUM) .AND. PRDZ) THEN U(NX,NY,1) = U(NX,NY,1) * + GEDGE0*( BD1(NY,1) + BD4(NX,1) ) * + GEDGE1*( BD1(NYM1,1) + BD4(NXM1,1) ) * + GEDGE2*( BD1(NYM2,1) + BD4(NXM2,1) ) * + GEDGE3*( BD1(NYM3,1) + BD4(NXM3,1) ) * + GEDGE4*( BD1(NY,2) + BD4(NX,2) * + BD1(NY,NZM1) + BD4(NX,NZM1) ) * + GEDGE5*( BD1(NYM1,2) + BD4(NXM1,2) * + BD1(NYM1,NZM1) + BD4(NXM1,NZM1) ) ENDIF C C C ---- C EXIT C ---- C 999 CONTINUE RETURN END SUBROUTINE MDALG3 (A, B, C, D, BCTY, U, LDXU, LDYU, IL, IR, * JL, JR, KL, KR, WORK) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C C DECEMBER 1985 (REVISED APRIL 1987) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C MDALG3 IMPLEMENTS THE MATRIX DECOMPOSITION ALGORITHM (FOURIER C METHOD FOR A NINE-POINT DIFFERENCE OPERATOR ON A THREE-DIMENSIONAL C RECTANGULAR GRID. C C C P A R A M E T E R S C ------------------- C C A, B, C, D REAL SCALARS (INPUT) C GIVE THE BASIC FINITE DIFFERENCE STENCIL THAT IS C USED TO APPROXIMATE THE PDE. C C D C D C C B C C D C D C C C B C C B A B U = RIGHT HAND SIDE C C B C C C D C D C C B C C D C D C C BCTY INTEGER ARRAY OF SIZE 6 (INPUT) C SEE HFFT3A. C C U REAL ARRAY OF SIZE LDXU BY LDYU BY KR (INPUT/OUTPUT) C ON INPUT, U CONTAINS THE RIGHT HAND SIDE OF THE C DISCRETE APPROXIMATION TO THE PDE FOR EACH POINT C POINT AT WHICH THE SOLUTION IS TO BE DETERMINED, C I.E., (I,J,K), I=IL,..,IR, J=JL,..,JR, K=KL,..,KR. C ON OUTPUT, THESE VALUES ARE REPLACED BY THE COMPUTED C SOLUTION. C C LDXU INTEGER SCALAR (INPUT) C THE LEADING DIMENSION OF THE ARRAY U EXACTLY AS C SPECIFIED IN THE CALLING PROGRAM. C C LDYU INTEGER SCALAR (INPUT) C THE SECOND DIMENSION OF THE ARRAY U EXACTLY AS C SPECIFIED IN THE CALLING PROGRAM. C C IL, IR, INTEGER SCALARS (INPUT) C JL, JR, GIVES THE SUBSET OF GRID POINTS AT WHICH THE C KL, KR SOLUTION IS TO BE DETERMINED, I.E., THAT SET C OF INDICES (I,J) WITH I=IL,..,IR, J=JL,..,JR, AND C K=KL,..,KR. C C WORK REAL ARRAY OF SIZE N*L + 5N + 5M + 3L + N/2 + L/2 +30 C HERE N=IR-IL+1, M=JR-JL+1, AND L=KR-KL+1. THE LENGTH C OF THIS ARRAY MAY BE REDUCED BY M WHEN THE SOLUTION C IS NOT PERIODIC IN Y, AND MAY BE REDUCED BY 4M WHEN C THE COEFFICIENT OF U IN THE PDE IS .LE. 0. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C INTEGER BCTY(6), LDXU, LDYU, IL, IR, JL, JR, KL, KR REAL * A, B, C, D, U(LDXU,LDYU,*), WORK(*) C C ... LOCAL VARIABLES C INTEGER N, M, L, LOCWSA, LOCMWK, TOTAL C C ... LOCAL CONSTANTS C INTEGER DRCH, NEUM, PRDC, LEFT, RIGHT, TOP, BOTTOM, FRONT, BACK PARAMETER (DRCH=1, NEUM =2, PRDC=3, * LEFT=3, RIGHT=1, TOP =4, BOTTOM=2, FRONT=5, BACK=6) C C C -------------- C INITIALIZATION C -------------- C N = IR - IL + 1 M = JR - JL + 1 L = KR - KL + 1 LOC WSA = 1 + N*L LOC MWK = LOC WSA + (3*L + L/2 + 15) TOTAL = LOC MWK + 5*N + 5*M + N/2 + 14 C CALL FFTI(BCTY(BACK),BCTY(FRONT),L,WORK(LOCWSA)) C C C ------------------ C FORWARD TRANSFORMS C ------------------ C DO 100 J=JL,JR KK = 0 DO 25 I=IL,IR DO 25 K=KL,KR KK = KK + 1 WORK(KK) = U(I,J,K) 25 CONTINUE CALL FFTF(BCTY(BACK),BCTY(FRONT),WORK,L,1,L,1,N,WORK(LOCWSA)) KK = 0 DO 75 I=IL,IR DO 75 K=KL,KR KK = KK + 1 U(I,J,K) = WORK(KK) 75 CONTINUE 100 CONTINUE C C C -------------------------------- C SOLVE N TWO-DIMENSIONAL PROBLEMS C -------------------------------- C CALL EVDISC(BCTY(BACK),BCTY(FRONT),WORK(KL),L) DO 200 K=KL,KR EVA = A + B*WORK(K) EVB = B + C*WORK(K) EVC = C + D*WORK(K) CALL MDALG2(EVA,EVB,EVC,BCTY,U(1,1,K),LDXU,IL,IR,JL,JR, * WORK(LOCMWK)) 200 CONTINUE C C ------------------- C BACKWARD TRANSFORMS C ------------------- C DO 300 J=JL,JR KK = 0 DO 225 I=IL,IR DO 225 K=KL,KR KK = KK + 1 WORK(KK) = U(I,J,K) 225 CONTINUE CALL FFTB(BCTY(BACK),BCTY(FRONT),WORK,L,1,L,1,N,WORK(LOCWSA)) KK = 0 DO 275 I=IL,IR DO 275 K=KL,KR KK = KK + 1 U(I,J,K) = WORK(KK) 275 CONTINUE 300 CONTINUE C C C ---- C EXIT C ---- C RETURN END SUBROUTINE REFL3 (KDIST, NX, NY, NZ, PRDX, PRDY, PRDZ, G, * LMXG, LMYG) C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C --------------- 4TH ORDER ACCURATE FAST DIRECT SOLUTION C PACKAGE : HFFT OF THE HELMHOLTZ EQUATION ON RECTANGULAR C --------------- DOMAINS IN TWO AND THREE DIMENSIONS C C INTERNAL MODULE C C RONALD F. BOISVERT C NATIONAL BUREAU OF STANDARDS C DECEMBER 1985 C C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C C REFL3 EXTENDS A THREE-DIMENSIONAL GRID FUNCTION SO THAT IT IS C DEFINED ONE GRID LINE OUTSIDE ITS ORIGINAL DOMAIN. THE EXTENSION C IS DONE BY REFLECTION THROUGH THE BOUNDARY EXCEPT WHERE C PERIODICITY IS SPECIFIED. C C C P A R A M E T E R S C ------------------- C C KDIST INTEGER SCALAR (INPUT) C INDICATES HOW THE GRID FUNCTION IS DEFINED. C POSSIBLE VALUES ARE C C 1 == FUNCTION DEFINED AT CENTER OF GRID SQUARES C 2 == FUNCTION DEFINED AT GRID POINTS C C NX INTEGER SCALAR (INPUT) C THE NUMBER OF GRID FUNCTION VALUES IN THE X DIRECTION. C C NY INTEGER SCALAR (INPUT) C THE NUMBER OF GRID FUNCTION VALUES IN THE Y DIRECTION. C C NZ INTEGER SCALAR (INPUT) C THE NUMBER OF GRID FUNCTION VALUES IN THE Z DIRECTION. C C PRDX LOGICAL SCALAR (INPUT) C .TRUE. IF THE SOLUTION IS TO BE EXTENDED PERIODICALLY C IN X. C C PRDY LOGICAL SCALAR (INPUT) C .TRUE. IF THE SOLUTION IS TO BE EXTENDED PERIODICALLY C IN Y. C C PRDZ LOGICAL SCALAR (INPUT) C .TRUE. IF THE SOLUTION IS TO BE EXTENDED PERIODICALLY C IN Z. C C G REAL ARRAY OF SIZE LMXG+1 BY LMYG+1 BY NZ+1 (INPUT/OUTPUT) C ON INPUT, THE GRID FUNCTION OCCUPIES G(I,J,K), I=1,..,NX, C J=1,..,NY, K=1,..,NZ. C ON OUTPUT, THE FUNCTION HAS BEEN EXTENDED TO INCLUDE THE C POINTS G(0,J,K), G(NX+1,J,K), G(I,0,K), G(I,NY+1,K), C G(I,J,0), G(I,J,NZ+1), I=0,..,NX+1, J=0,..,NY+1, C K=0,..,NZ+1. C C LMXG INTEGER SCALAR (INPUT) C THE UPPER LIMIT OF THE FIRST DIMENSION OF THE ARRAY G. C MUST BE SET TO LDXU-1, WHERE LDXU IS THE ACTUAL LENGTH OF C THE FIRST DIMENSION OF G AS DECLARED IN THE CALLING C PROGRAM. C C LMYG INTEGER SCALAR (INPUT) C THE UPPER LIMIT OF THE SECOND DIMENSION OF THE ARRAY G. C MUST BE SET TO LDYU-1, WHERE LDYU IS THE ACTUAL LENGTH OF C THE SECOND DIMENSION OF G AS DECLARED IN THE CALLING C PROGRAM. C C C ------------ C DECLARATIONS C ------------ C C ... PARAMETERS C LOGICAL PRDX, PRDY, PRDZ INTEGER KDIST, NX, NY, NZ, LMXG, LMYG REAL * G(0:LMXG,0:LMYG,0:*) C C C --------------- C INITIAL