SUBROUTINE PWSCRT (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC, POI00001 1 BDD,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE PWSCRT SOLVES THE STANDARD FIVE-POINT FINITE C DIFFERENCE APPROXIMATION TO THE HELMHOLTZ EQUATION IN CARTESIAN C COORDINATES: C C (D/DX)(D/DX)U + (D/DY)(D/DY)U + LAMBDA*U = F(X,Y). C C THE ARGUMENTS ARE DEFINED AS: C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C INTL C = 0 ON INITIAL ENTRY TO PWSCRT OR IF N OR NBDCND ARE CHANGED C FROM PREVIOUS CALL. C = 1 IF N AND NBDCND ARE UNCHANGED FROM PREVIOUS CALL TO PWSCRT. C C NOTE: A CALL WITH INTL = 1 IS ABOUT 1 PER CENT FASTER THAN A C CALL WITH INTL = 0 . C C A,B C THE RANGE OF X, I.E., A .LE. X .LE. B. A MUST BE LESS THAN B. C C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (A,B) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C X-DIRECTION GIVEN BY X(I) = A+(I-1)DX FOR I = 1,2,...,M+1, C WHERE DX = (B-A)/M IS THE PANEL WIDTH. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT X = A AND X = B. C C = 0 IF THE SOLUTION IS PERIODIC IN X, I.E., U(1,J) = U(M+1,J). C = 1 IF THE SOLUTION IS SPECIFIED AT X = A AND X = B. C = 2 IF THE SOLUTION IS SPECIFIED AT X = A AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO X IS SPECIFIED AT X = B. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = A AND X = B. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X IS C SPECIFIED AT X = A AND THE SOLUTION IS SPECIFIED AT X = B. C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X AT X = A. C WHEN MBDCND = 3 OR 4, C C BDA(J) = (D/DX)U(A,Y(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A DUMMY VARIABLE. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO X AT X = B. C WHEN MBDCND = 2 OR 3, C C BDB(J) = (D/DX)U(B,Y(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE BDB IS A DUMMY VARIABLE. C C C,D C THE RANGE OF Y, I.E., C .LE. Y .LE. D. C MUST BE LESS THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (C,D) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C Y-DIRECTION GIVEN BY Y(J) = C+(J-1)DY FOR J = 1,2,...,N+1, WHERE C DY = (D-C)/N IS THE PANEL WIDTH. N MUST BE OF THE FORM C (2**P)(3**Q)(5**R) WHERE P, Q, AND R ARE ANY NON-NEGATIVE C INTEGERS. N MUST BE GREATER THAN 2. C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT Y = C AND Y = D. C C = 0 IF THE SOLUTION IS PERIODIC IN Y, I.E., U(I,1) = U(I,N+1). C = 1 IF THE SOLUTION IS SPECIFIED AT Y = C AND Y = D. C = 2 IF THE SOLUTION IS SPECIFIED AT Y = C AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO Y IS SPECIFIED AT Y = D. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = C AND Y = D. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y IS C SPECIFIED AT Y = C AND THE SOLUTION IS SPECIFIED AT Y = D. C C BDC C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y AT Y = C. C WHEN NBDCND = 3 OR 4, C C BDC(I) = (D/DY)U(X(I),C), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS A DUMMY VARIABLE. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Y AT Y = D. C WHEN NBDCND = 2 OR 3, C C BDD(I) = (D/DY)U(X(I),D), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, PWSCRT WILL C ATTEMPT TO FIND A SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY WHICH SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY). C FOR I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(X(I),Y(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ --------- -------- C C 0 F(A,Y(J)) F(A,Y(J)) C 1 U(A,Y(J)) U(B,Y(J)) C 2 U(A,Y(J)) F(B,Y(J)) J = 1,2,...,N+1 C 3 F(A,Y(J)) F(B,Y(J)) C 4 F(A,Y(J)) U(B,Y(J)) C C C NBDCND F(I,1) F(I,N+1) C ------ --------- -------- C C 0 F(X(I),C) F(X(I),C) C 1 U(X(I),C) U(X(I),D) C 2 U(X(I),C) F(X(I),D) I = 1,2,...,M+1 C 3 F(X(I),C) F(X(I),D) C 4 F(X(I),C) U(X(I),D) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING PWSCRT. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 6(N+1)+8(M+1). C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (X(I),Y(J)), I = 1,2,...,M+1, C J = 1,2,...,N+1 . C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE BOUNDARY CONDITIONS C IS SPECIFIED FOR A POISSON EQUATION (LAMBDA = 0), A SOLUTION MAY C NOT EXIST. PERTRB IS A CONSTANT, CALCULATED AND SUBTRACTED FROM C F, WHICH ENSURES THAT A SOLUTION EXISTS. PWSCRT THEN COMPUTES C THIS SOLUTION, WHICH IS A LEAST SQUARES SOLUTION TO THE ORIGINAL C APPROXIMATION. THIS SOLUTION IS NOT UNIQUE AND IS UNNORMALIZED. C THE VALUE OF PERTRB SHOULD BE SMALL COMPARED TO THE RIGHT SIDE C F. OTHERWISE , A SOLUTION IS OBTAINED TO AN ESSENTIALLY C DIFFERENT PROBLEM. THIS COMPARISON SHOULD ALWAYS BE MADE TO C INSURE THAT A MEANINGFUL SOLUTION HAS BEEN OBTAINED. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 6, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 A .GE. B. C = 2 MBDCND .LT. 0 OR MBDCND .GT. 4 . C = 3 C .GE. D. C = 4 N .NE. (2**P)(3**Q)(5**R) OR N .LE. 2 . C = 5 NBDCND .LT. 0 OR NBDCND .GT. 4 . C = 6 LAMBDA .GT. 0 . C = 7 IDIMF .LT. M+1 . C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO PWSCRT, THE USER SHOULD TEST IERROR AFTER THE CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C PWSCRT WILL BE CALLED AGAIN WITH INTL = 1 . C C C DIMENSION F(IDIMF,1) DIMENSION BDA(1) ,BDB(1) ,BDC(1) ,BDD(1) , 1 W(1) C C CHECK FOR INVALID PARAMETERS. C IERROR = 0 IF (A .GE. B) IERROR = 1 IF (MBDCND.LT.0 .OR. MBDCND.GT.4) IERROR = 2 IF (C .GE. D) IERROR = 3 IF (NCHECK(N).GT.1 .OR. N.LE.2) IERROR = 4 IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 5 IF (IDIMF .LT. M+1) IERROR = 7 IF (IERROR .NE. 0) RETURN NPEROD = NBDCND NBY2 = 0 DELTAX = (B-A)/FLOAT(M) TWDELX = 2./DELTAX DELXSQ = 1./DELTAX**2 DELTAY = (D-C)/FLOAT(N) TWDELY = 2./DELTAY DELYSQ = 1./DELTAY**2 NP = NBDCND+1 NP1 = N+1 MP = MBDCND+1 MP1 = M+1 NSTART = 1 NSTOP = N NSKIP = 1 GO TO (104,101,102,103,104),NP 101 NSTART = 2 GO TO 104 102 NSTART = 2 103 NSTOP = NP1 NSKIP = 2 104 NUNK = NSTOP-NSTART+1 C C ENTER BOUNDARY DATA FOR X-BOUNDARIES. C MSTART = 1 MSTOP = M MSKIP = 1 GO TO (117,105,106,109,110),MP 105 MSTART = 2 GO TO 107 106 MSTART = 2 MSTOP = MP1 MSKIP = 2 107 DO 108 J=NSTART,NSTOP F(2,J) = F(2,J)-F(1,J)*DELXSQ 108 CONTINUE GO TO 112 109 MSTOP = MP1 MSKIP = 2 110 DO 111 J=NSTART,NSTOP F(1,J) = F(1,J)+BDA(J)*TWDELX 111 CONTINUE 112 GO TO (113,115),MSKIP 113 DO 114 J=NSTART,NSTOP F(M,J) = F(M,J)-F(MP1,J)*DELXSQ 114 CONTINUE GO TO 117 115 DO 116 J=NSTART,NSTOP F(MP1,J) = F(MP1,J)-BDB(J)*TWDELX 116 CONTINUE 117 MUNK = MSTOP-MSTART+1 C C ENTER BOUNDARY DATA FOR Y-BOUNDARIES. C GO TO (127,118,118,120,120),NP 118 DO 119 I=MSTART,MSTOP F(I,2) = F(I,2)-F(I,1)*DELYSQ 119 CONTINUE GO TO 122 120 DO 121 I=MSTART,MSTOP F(I,1) = F(I,1)+BDC(I)*TWDELY 121 CONTINUE 122 GO TO (123,125),NSKIP 123 DO 124 I=MSTART,MSTOP F(I,N) = F(I,N)-F(I,NP1)*DELYSQ 124 CONTINUE GO TO 127 125 DO 126 I=MSTART,MSTOP F(I,NP1) = F(I,NP1)-BDD(I)*TWDELY 126 CONTINUE C C MULTIPLY RIGHT SIDE BY DELTAY**2. C 127 DELYSQ = DELTAY*DELTAY DO 129 I=MSTART,MSTOP DO 128 J=NSTART,NSTOP F(I,J) = F(I,J)*DELYSQ 128 CONTINUE 129 CONTINUE C C DEFINE THE A,B,C COEFFICIENTS IN W-ARRAY. C ID1 = 6*NUNK+5*MUNK ID2 = ID1+MUNK ID3 = ID2+MUNK S = DELYSQ*DELXSQ ST2 = 2.*S DO 130 I=1,MUNK J = ID1+I W(J) = S J = ID2+I W(J) = -ST2+ELMBDA*DELYSQ J = ID3+I W(J) = S 130 CONTINUE GO TO (134,134,131,132,133),MP 131 W(ID2) = ST2 GO TO 134 132 W(ID2) = ST2 133 W(ID3+1) = ST2 134 CONTINUE IF (NBDCND .NE. 4) GO TO 138 C C REVERSE COLUMNS TO CHANGE DERIVATIVE SPECIFIED-SOLUTION SPECIFIED C BOUNDARY CONDITIONS TO SOLUTION SPECIFIED-DERIVATIVE SPECIFIED. C IREV = 1 NPEROD = 2 NBY2 = N/2 135 DO 137 J=1,NBY2 MSKIP = N+1-J DO 136 I=MSTART,MSTOP A1 = F(I,J) F(I,J) = F(I,MSKIP) F(I,MSKIP) = A1 136 CONTINUE 137 CONTINUE GO TO (147,150),IREV 138 CONTINUE PERTRB = 0. IF (ELMBDA) 147,140,139 139 IERROR = 6 GO TO 147 140 IF ((NBDCND.EQ.0 .OR. NBDCND.EQ.3) .AND. 1 (MBDCND.EQ.0 .OR. MBDCND.EQ.3)) GO TO 141 GO TO 147 C C FOR SINGULAR PROBLEMS MUST ADJUST DATA TO INSURE THAT A SOLUTION C WILL EXIST. C 141 A1 = 1. A2 = 1. IF (NBDCND .EQ. 3) A2 = 2. IF (MBDCND .EQ. 3) A1 = 2. S1 = 0. MSP1 = MSTART+1 MSTM1 = MSTOP-1 NSP1 = NSTART+1 NSTM1 = NSTOP-1 DO 143 J=NSP1,NSTM1 S = 0. DO 142 I=MSP1,MSTM1 S = S+F(I,J) 142 CONTINUE S1 = S1+S*A1+F(MSTART,J)+F(MSTOP,J) 143 CONTINUE S1 = A2*S1 S = 0. DO 144 I=MSP1,MSTM1 S = S+F(I,NSTART)+F(I,NSTOP) 144 CONTINUE S1 = S1+S*A1+F(MSTART,NSTART)+F(MSTART,NSTOP)+F(MSTOP,NSTART)+ 1 F(MSTOP,NSTOP) S = (2.+FLOAT(NUNK-2)*A2)*(2.+FLOAT(MUNK-2)*A1) PERTRB = S1/S DO 146 J=NSTART,NSTOP DO 145 I=MSTART,MSTOP F(I,J) = F(I,J)-PERTRB 145 CONTINUE 146 CONTINUE PERTRB = PERTRB/DELYSQ C C SOLVE THE EQUATION. C 147 CALL POIS (INTL,NPEROD,NUNK,MBDCND,MUNK,W(ID1+1),W(ID2+1), 1 W(ID3+1),IDIMF,F(MSTART,NSTART),IERR1,W) IF (NBY2 .EQ. 0) GO TO 148 IREV = 2 GO TO 135 148 CONTINUE C C FILL IN IDENTICAL VALUES WHEN HAVE PERIODIC BOUNDARY CONDITIONS. C IF (NBDCND .NE. 0) GO TO 150 DO 149 I=MSTART,MSTOP F(I,NP1) = F(I,1) 149 CONTINUE 150 IF (MBDCND .NE. 0) GO TO 152 DO 151 J=NSTART,NSTOP F(MP1,J) = F(1,J) 151 CONTINUE IF (NBDCND .EQ. 0) F(MP1,NP1) = F(1,NP1) 152 CONTINUE RETURN END SUBROUTINE PWSPLR (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC, POI00418 1 BDD,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE PWSPLR SOLVES A FINITE DIFFERENCE APPROXIMATION TO THE C HELMHOLTZ EQUATION IN POLAR COORDINATES: C C (1/R)(D/DR)(R(D/DR)U) + (1/R**2)(D/DTHETA)(D/DTHETA)U C C + LAMBDA*U = F(R,THETA). C C THE ARGUMENTS ARE DEFINED AS : C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C INTL C = 0 ON INITIAL ENTRY TO PWSPLR OR IF N OR NBDCND ARE CHANGED C FROM PREVIOUS CALL. C = 1 IF N AND NBDCND ARE UNCHANGED FROM PREVIOUS CALL TO PWSPLR. C C NOTE: A CALL WITH INTL = 1 IS ABOUT 1 PER CENT FASTER THAN A C CALL WITH INTL = 0 . C C A,B C THE RANGE OF R, I.E., A .LE. R .LE. B. A MUST BE LESS THAN B C AND A MUST BE NON-NEGATIVE. C C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (A,B) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C R-DIRECTION GIVEN BY R(I) = A+(I-1)DR, FOR I = 1,2,...,M+1, C WHERE DR = (B-A)/M IS THE PANEL WIDTH. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION AT R = A AND R = B. C C = 1 IF THE SOLUTION IS SPECIFIED AT R = A AND R = B. C = 2 IF THE SOLUTION IS SPECIFIED AT R = A AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO R IS SPECIFIED AT R = B. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = A (SEE NOTE BELOW) AND R = B. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = A (SEE NOTE BELOW) AND THE SOLUTION IS C SPECIFIED AT R = B. C = 5 IF THE SOLUTION IS UNSPECIFIED AT R = A = 0 AND THE C SOLUTION IS SPECIFIED AT R = B. C = 6 IF THE SOLUTION IS UNSPECIFIED AT R = A = 0 AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS SPECIFIED C AT R = B. C C NOTE: IF A = 0, DO NOT USE MBDCND = 3 OR 4, BUT INSTEAD USE C MBDCND = 1,2,5, OR 6 . C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = A. C WHEN MBDCND = 3 OR 4, C C BDA(J) = (D/DR)U(A,THETA(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A DUMMY VARIABLE. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = B. C WHEN MBDCND = 2,3, OR 6, C C BDB(J) = (D/DR)U(B,THETA(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDB IS A DUMMY VARIABLE. C C C,D C THE RANGE OF THETA, I.E., C .LE. THETA .LE. D. C MUST BE LESS C THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (C,D) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C THETA-DIRECTION GIVEN BY THETA(J) = C+(J-1)DTHETA FOR C J = 1,2,...,N+1, WHERE DTHETA = (D-C)/N IS THE PANEL WIDTH. N C MUST BE IN THE FORM (2**P)(3**Q)(5**R) WHERE P, Q, AND R ARE ANY C NON-NEGATIVE INTEGERS. N MUST BE GREATER THAN 2. C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT THETA = C AND C AND THETA = D. C C = 0 IF THE SOLUTION IS PERIODIC IN THETA, I.E., C U(I,1) = U(I,N+1). C = 1 IF THE SOLUTION IS SPECIFIED AT THETA = C AND THETA = D C (SEE NOTE BELOW). C = 2 IF THE SOLUTION IS SPECIFIED AT THETA = C AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = D (SEE NOTE BELOW). C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = C AND THE SOLUTION IS SPECIFIED AT C THETA = D (SEE NOTE BELOW). C C NOTE: WHEN NBDCND = 1,2, OR 4, DO NOT USE MBDCND = 5 OR 6 C (THE FORMER INDICATES THAT THE SOLUTION IS SPECIFIED AT C R = 0, THE LATTER INDICATES THE SOLUTION IS UNSPECIFIED C AT R = 0). USE INSTEAD MBDCND = 1 OR 2 . C C BDC C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = C. WHEN NBDCND = 3 OR 4, C C BDC(I) = (D/DTHETA)U(R(I),C), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS A DUMMY VARIABLE. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = D. WHEN NBDCND = 2 OR 3, C C BDD(I) = (D/DTHETA)U(R(I),D), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, PWSPLR WILL C ATTEMPT TO FIND A SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY). C FOR I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(R(I),THETA(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ ------------- ------------- C C 1 U(A,THETA(J)) U(B,THETA(J)) C 2 U(A,THETA(J)) F(B,THETA(J)) C 3 F(A,THETA(J)) F(B,THETA(J)) C 4 F(A,THETA(J)) U(B,THETA(J)) J = 1,2,...,N+1 C 5 F(0,0) U(B,THETA(J)) C 6 F(0,0) F(B,THETA(J)) C C NBDCND F(I,1) F(I,N+1) C ------ --------- --------- C C 0 F(R(I),C) F(R(I),C) C 1 U(R(I),C) U(R(I),D) C 2 U(R(I),C) F(R(I),D) I = 1,2,...,M+1 C 3 F(R(I),C) F(R(I),D) C 4 F(R(I),C) U(R(I),D) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING PWSPLR. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 6(N+1)+8(M+1). C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (R(I),THETA(J)), C I = 1,2,...,M+1, J = 1,2,...,N+1 . C C PERTRB C IF A COMBINATION OF PERIODIC, DERIVATIVE, OR UNSPECIFIED C BOUNDARY CONDITIONS IS SPECIFIED FOR A POISSON EQUATION C (LAMBDA = 0), A SOLUTION MAY NOT EXIST. PERTRB IS A CONSTANT, C CALCULATED AND SUBTRACTED FROM F, WHICH ENSURES THAT A SOLUTION C EXISTS. PWSPLR THEN COMPUTES THIS SOLUTION, WHICH IS A LEAST C SQUARES SOLUTION TO THE ORIGINAL APPROXIMATION. THIS SOLUTION C IS NOT UNIQUE AND IS UNNORMALIZED. THE VALUE OF PERTRB SHOULD C BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE , A SOLUTION C IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. THIS COMPARISON C SHOULD ALWAYS BE MADE TO INSURE THAT A MEANINGFUL SOLUTION HAS C BEEN OBTAINED C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 11, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 A .LT. 0 . C = 2 A .GE. B. C = 3 MBDCND .LT. 1 OR MBDCND .GT. 6 . C = 4 C .GE. D. C = 5 N .NE. (2**P)(3**Q)(5**R) OR N .LE. 2 . C = 6 NBDCND .LT. 0 OR .GT. 4 . C = 7 A = 0, MBDCND = 3 OR 4 . C = 8 A .GT. 0, MBDCND .GE. 5 . C = 9 MBDCND .GE. 5, NBDCND .NE. 0 AND NBDCND .NE. 3 . C = 10 IDIMF .LT. M+1 . C = 11 LAMBDA .GT. 0 . C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO PWSPLR, THE USER SHOULD TEST IERROR AFTER THE CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C PWSPLR WILL BE CALLED AGAIN WITH INTL = 1 . C C C DIMENSION F(IDIMF,1) DIMENSION BDA(1) ,BDB(1) ,BDC(1) ,BDD(1) , 1 W(1) C C CHECK FOR INVALID PARAMETERS. C IERROR = 0 IF (A .LT. 0.) IERROR = 1 IF (A .GE. B) IERROR = 2 IF (MBDCND.LE.0 .OR. MBDCND.GE.7) IERROR = 3 IF (C .GE. D) IERROR = 4 IF (NCHECK(N).GT.1 .OR. N.LE.2) IERROR = 5 IF (NBDCND.LE.-1 .OR. NBDCND.GE.5) IERROR = 6 IF (A.EQ.0. .AND. (MBDCND.EQ.3 .OR. MBDCND.EQ.4)) IERROR = 7 IF (A.GT.0. .AND. MBDCND.GE.5) IERROR = 8 IF (MBDCND.GE.5 .AND. NBDCND.NE.0 .AND. NBDCND.NE.3) IERROR = 9 IF (IDIMF .LT. M+1) IERROR = 10 IF (IERROR .NE. 0) RETURN MP1 = M+1 DELTAR = (B-A)/FLOAT(M) DLRBY2 = DELTAR/2. DLRSQ = DELTAR**2 NP1 = N+1 DELTHT = (D-C)/FLOAT(N) DLTHSQ = DELTHT**2 NP = NBDCND+1 C C DEFINE RANGE OF INDICES I AND J FOR UNKNOWNS U(I,J). C MSTART = 2 MSTOP = MP1 GO TO (101,105,102,103,104,105),MBDCND 101 MSTOP = M GO TO 105 102 MSTART = 1 GO TO 105 103 MSTART = 1 104 MSTOP = M 105 MUNK = MSTOP-MSTART+1 NSTART = 1 NSTOP = N GO TO (109,106,107,108,109),NP 106 NSTART = 2 GO TO 109 107 NSTART = 2 108 NSTOP = NP1 109 NUNK = NSTOP-NSTART+1 C C DEFINE A,B,C COEFFICIENTS IN W-ARRAY. C ID1 = 6*NUNK+5*MUNK ID2 = ID1+MUNK ID3 = ID2+MUNK ID4 = ID3+MUNK ID5 = ID1-5*MUNK ID6 = ID5+MUNK A1 = 2./DLRSQ IJ = 0 IF (MBDCND.EQ.3 .OR. MBDCND.EQ.4) IJ = 1 DO 110 I=1,MUNK R = A+FLOAT(I-IJ)*DELTAR J = ID5+I W(J) = R J = ID6+I W(J) = 1./R**2 J = ID1+I W(J) = (R-DLRBY2)/(R*DLRSQ) J = ID3+I W(J) = (R+DLRBY2)/(R*DLRSQ) J = ID2+I W(J) = -A1+ELMBDA 110 CONTINUE GO TO (114,111,112,113,114,111),MBDCND 111 W(ID2) = A1 GO TO 114 112 W(ID2) = A1 113 W(ID3+1) = A1 114 CONTINUE C C ENTER BOUNDARY DATA FOR R-BOUNDARIES. C GO TO (115,115,117,117,119,119),MBDCND 115 A1 = W(ID1+1) DO 116 J=NSTART,NSTOP F(2,J) = F(2,J)-A1*F(1,J) 116 CONTINUE GO TO 119 117 A1 = 2.*DELTAR*W(ID1+1) DO 118 J=NSTART,NSTOP F(1,J) = F(1,J)+A1*BDA(J) 118 CONTINUE 119 GO TO (120,122,122,120,120,122),MBDCND 120 A1 = W(ID4) DO 121 J=NSTART,NSTOP F(M,J) = F(M,J)-A1*F(MP1,J) 121 CONTINUE GO TO 124 122 A1 = 2.*DELTAR*W(ID4) DO 123 J=NSTART,NSTOP F(MP1,J) = F(MP1,J)-A1*BDB(J) 123 CONTINUE C C ENTER BOUNDARY DATA FOR THETA-BOUNDARIES. C 124 A1 = 1./DLTHSQ L = ID5-MSTART+1 LP = ID6-MSTART+1 GO TO (134,125,125,127,127),NP 125 DO 126 I=MSTART,MSTOP J = I+LP F(I,2) = F(I,2)-A1*W(J)*F(I,1) 126 CONTINUE GO TO 129 127 A1 = 2./DELTHT DO 128 I=MSTART,MSTOP J = I+LP F(I,1) = F(I,1)+A1*W(J)*BDC(I) 128 CONTINUE 129 A1 = 1./DLTHSQ GO TO (134,130,132,132,130),NP 130 DO 131 I=MSTART,MSTOP J = I+LP F(I,N) = F(I,N)-A1*W(J)*F(I,NP1) 131 CONTINUE GO TO 134 132 A1 = 2./DELTHT DO 133 I=MSTART,MSTOP J = I+LP F(I,NP1) = F(I,NP1)-A1*W(J)*BDD(I) 133 CONTINUE 134 CONTINUE C C ADJUST RIGHT SIDE OF EQUATION FOR UNKNOWN AT POLE WHEN HAVE C DERIVATIVE SPECIFIED BOUNDARY CONDITIONS. C IF (MBDCND.GE.5 .AND. NBDCND.EQ.3) 1 F(1,1) = F(1,1)-(BDD(2)-BDC(2))*4./(FLOAT(N)*DELTHT*DLRSQ) C C ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A C SOLUTION. C PERTRB = 0. IF (ELMBDA) 144,136,135 135 IERROR = 11 GO TO 144 136 IF (NBDCND.NE.0 .AND. NBDCND.NE.3) GO TO 144 S2 = 0. GO TO (144,144,137,144,144,138),MBDCND 137 W(ID5+1) = .5*(W(ID5+2)-DLRBY2) S2 = .25*DELTAR 138 A2 = 2. IF (NBDCND .EQ. 0) A2 = 1. J = ID5+MUNK W(J) = .5*(W(J-1)+DLRBY2) S = 0. DO 140 I=MSTART,MSTOP S1 = 0. IJ = NSTART+1 K = NSTOP-1 DO 139 J=IJ,K S1 = S1+F(I,J) 139 CONTINUE J = I+L S = S+(A2*S1+F(I,NSTART)+F(I,NSTOP))*W(J) 140 CONTINUE S2 = FLOAT(M)*A+DELTAR*(FLOAT((M-1)*(M+1))*.5+.25)+S2 S1 = (2.+A2*FLOAT(NUNK-2))*S2 IF (MBDCND .EQ. 3) GO TO 141 S2 = FLOAT(N)*A2*DELTAR/8. S = S+F(1,1)*S2 S1 = S1+S2 141 CONTINUE PERTRB = S/S1 DO 143 I=MSTART,MSTOP DO 142 J=NSTART,NSTOP F(I,J) = F(I,J)-PERTRB 142 CONTINUE 143 CONTINUE 144 CONTINUE C C MULTIPLY I-TH EQUATION THROUGH BY (R(I)*DELTHT)**2. C DO 146 I=MSTART,MSTOP K = I-MSTART+1 J = I+LP A1 = DLTHSQ/W(J) J = ID1+K W(J) = A1*W(J) J = ID2+K W(J) = A1*W(J) J = ID3+K W(J) = A1*W(J) DO 145 J=NSTART,NSTOP F(I,J) = A1*F(I,J) 145 CONTINUE 146 CONTINUE LP = NBDCND IF (LP .NE. 4) GO TO 150 C C REVERSE COLUMNS SO THE DERIVATIVE SPECIFIED-SOLUTION SPECIFIED C BOUNDARY CONDITIONS BECOME SOLUTION SPECIFIED-DERIVATIVE SPECIFIED C LP = 2 ID5 = 1 K = N/2 147 DO 149 J=1,K L = N+1-J DO 148 I=MSTART,MSTOP S = F(I,L) F(I,L) = F(I,J) F(I,J) = S 148 CONTINUE 149 CONTINUE GO TO (150,162),ID5 C C CALL POIS TO SOLVE THE SYSTEM OF EQUATIONS. C 150 CALL POIS (INTL,LP,NUNK,MBDCND,MUNK,W(ID1+1),W(ID2+1),W(ID3+1), 1 IDIMF,F(MSTART,NSTART),IERR1,W) IF (NBDCND .NE. 4) GO TO 151 ID5 = 2 GO TO 147 151 GO TO (162,162,162,162,153,152),MBDCND C C ADJUST THE SOLUTION AS NECESSARY FOR THE PROBLEMS WHERE A = 0. C 152 IF (ELMBDA .NE. 0.) GO TO 153 YPOLE = 0. GO TO 160 153 CONTINUE J = ID5+MUNK W(J) = W(ID2)/W(ID3) DO 154 IP=3,MUNK I = MUNK-IP+2 J = ID5+I IJ = ID1+I LP = ID2+I K = ID3+I W(J) = W(IJ)/(W(LP)-W(K)*W(J+1)) 154 CONTINUE W(ID5+1) = -.5*DLTHSQ/(W(ID2+1)-W(ID3+1)*W(ID5+2)) DO 155 I=2,MUNK J = ID5+I W(J) = -W(J)*W(J-1) 155 CONTINUE S = 0. DO 156 J=NSTART,NSTOP S = S+F(2,J) 156 CONTINUE A2 = NUNK IF (NBDCND .EQ. 0) GO TO 157 S = S-.5*(F(2,NSTART)+F(2,NSTOP)) A2 = A2-1. 157 YPOLE = (.25*DLRSQ*F(1,1)-S/A2)/(W(ID5+1)-1.+ELMBDA*DLRSQ*.25) DO 159 I=MSTART,MSTOP K = L+I DO 158 J=NSTART,NSTOP F(I,J) = F(I,J)+YPOLE*W(K) 158 CONTINUE 159 CONTINUE 160 DO 161 J=1,NP1 F(1,J) = YPOLE 161 CONTINUE 162 CONTINUE IF (NBDCND .NE. 0) GO TO 164 DO 163 I=MSTART,MSTOP F(I,NP1) = F(I,1) 163 CONTINUE 164 CONTINUE RETURN END SUBROUTINE PWSCYL (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC, POI00931 1 BDD,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE PWSCYL SOLVES A FINITE DIFFERENCE APPROXIMATION TO THE C MODIFIED HELMHOLTZ EQUATION IN CYLINDRICAL COORDINATES C C (1/R)(D/DR)(R(D/DR)U) + (D/DZ)(D/DZ)U C C + (LAMBDA/R**2)U = F(R,Z) C C THIS TWO DIMENSIONAL MODIFIED HELMHOLTZ EQUATION RESULTS FROM C THE FOURIER TRANSFORM OF THE THREE DIMENSIONAL POISSON EQUATION C C C THE ARGUMENTS ARE DEFINED AS: C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C INTL C = 0 ON INITIAL ENTRY TO PWSCYL OR IF N AND NBDCND ARE CHANGED C FROM PREVIOUS CALL. C = 1 IF N AND NBDCND ARE UNCHANGED FROM PREVIOUS CALL TO PWSCYL. C C NOTE: A CALL WITH INTL = 1 IS ABOUT 1 PERCENT FASTER THAN A C CALL WITH INTL = 0 . C C A,B C THE RANGE OF R, I.E., A .LE. R .LE. B. A MUST BE LESS THAN B C AND A MUST BE NON-NEGATIVE. C C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (A,B) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C R-DIRECTION GIVEN BY R(I) = A+(I-1)DR, FOR I = 1,2,...,M+1, C WHERE DR = (B-A)/M IS THE PANEL WIDTH. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT R = A AND R = B. C C = 1 IF THE SOLUTION IS SPECIFIED AT R = A AND R = B. C = 2 IF THE SOLUTION IS SPECIFIED AT R = A AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO R IS SPECIFIED AT R = B. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = A (SEE NOTE BELOW) AND R = B. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = A (SEE NOTE BELOW) AND THE SOLUTION IS C SPECIFIED AT R = B. C = 5 IF THE SOLUTION IS UNSPECIFIED AT R = A = 0 AND THE C SOLUTION IS SPECIFIED AT R = B. C = 6 IF THE SOLUTION IS UNSPECIFIED AT R = A = 0 AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS SPECIFIED C AT R = B. C C NOTE: IF A = 0, DO NOT USE MBDCND = 3 OR 4, BUT INSTEAD USE C MBDCND = 1,2,5, OR 6 . C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = A. C WHEN MBDCND = 3 OR 4, C C BDA(J) = (D/DR)U(A,Z(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A DUMMY VARIABLE. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = B. C WHEN MBDCND = 2,3, OR 6, C C BDB(J) = (D/DR)U(B,Z(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDB IS A DUMMY VARIABLE. C C C,D C THE RANGE OF Z, I.E., C .LE. Z .LE. D. C MUST BE LESS THAN D. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (C,D) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C Z-DIRECTION GIVEN BY Z(J) = C+(J-1)DZ, FOR J = 1,2,...,N+1, C WHERE DZ = (D-C)/N IS THE PANEL WIDTH. N MUST BE OF THE FORM C (2**P)(3**Q)(5**R) WHERE P, Q, AND R ARE ANY NON-NEGATIVE C INTEGERS. N MUST BE GREATER THAN 2 . C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS AT Z = C AND Z = D. C C = 0 IF THE SOLUTION IS PERIODIC IN Z, I.E., U(I,1) = U(I,N+1). C = 1 IF THE SOLUTION IS SPECIFIED AT Z = C AND Z = D. C = 2 IF THE SOLUTION IS SPECIFIED AT Z = C AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO Z IS SPECIFIED AT Z = D. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z IS C SPECIFIED AT Z = C AND Z = D. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z IS C SPECIFIED AT Z = C AND THE SOLUTION IS SPECIFIED AT Z = D. C C BDC C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z AT Z = C. C WHEN NBDCND = 3 OR 4, C C BDC(I) = (D/DZ)U(R(I),C), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS A DUMMY VARIABLE. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO Z AT Z = D. C WHEN NBDCND = 2 OR 3, C C BDD(I) = (D/DZ)U(R(I),D), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, PWSCYL WILL C ATTEMPT TO FIND A SOLUTION. LAMBDA MUST BE ZERO WHEN C MBDCND = 5 OR 6 . C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY DATA (IF ANY). FOR C I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(R(I),Z(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ --------- --------- C C 1 U(A,Z(J)) U(B,Z(J)) C 2 U(A,Z(J)) F(B,Z(J)) C 3 F(A,Z(J)) F(B,Z(J)) J = 1,2,...,N+1 C 4 F(A,Z(J)) U(B,Z(J)) C 5 F(0,Z(J)) U(B,Z(J)) C 6 F(0,Z(J)) F(B,Z(J)) C C NBDCND F(I,1) F(I,N+1) C ------ --------- --------- C C 0 F(R(I),C) F(R(I),C) C 1 U(R(I),C) U(R(I),D) C 2 U(R(I),C) F(R(I),D) I = 1,2,...,M+1 C 3 F(R(I),C) F(R(I),D) C 4 F(R(I),C) U(R(I),D) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING PWSCYL. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 6(N+1)+8(M+1). C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (R(I),Z(J)), I = 1,2,...,M+1, C J = 1,2,...,N+1 . C C PERTRB C IF ONE SPECIFIES A COMBINATION OF PERIODIC, DERIVATIVE, AND C UNSPECIFIED BOUNDARY CONDITIONS FOR A POISSON EQUATION C (LAMBDA = 0), A SOLUTION MAY NOT EXIST. PERTRB IS A CONSTANT, C CALCULATED AND SUBTRACTED FROM F, WHICH ENSURES THAT A SOLUTION C EXISTS. PWSCYL THEN COMPUTES THIS SOLUTION, WHICH IS A LEAST C SQUARES SOLUTION TO THE ORIGINAL APPROXIMATION. THIS SOLUTION C IS NOT UNIQUE AND IS UNNORMALIZED. THE VALUE OF PERTRB SHOULD C BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE , A SOLUTION C IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. THIS COMPARISON C SHOULD ALWAYS BE MADE TO INSURE THAT A MEANINGFUL SOLUTION HAS C BEEN OBTAINED C C IERROR C AN ERROR FLAG WHICH INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 11, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 A .LT. 0 . C = 2 A .GE. B. C = 3 MBDCND .LT. 1 OR MBDCND .GT. 6 . C = 4 C .GE. D. C = 5 N .NE. (2**P)(3**Q)(5**R) OR N .LE. 2 . C = 6 NBDCND .LT. 0 OR NBDCND .GT. 4 . C = 7 A = 0, MBDCND = 3 OR 4 . C = 8 A .GT. 0, MBDCND .GE. 5 . C = 9 A = 0, LAMBDA .NE. 0, MBDCND .GE. 5 . C = 10 IDIMF .LT. M+1 . C = 11 LAMBDA .GT. 0 . C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO PWSCYL, THE USER SHOULD TEST IERROR AFTER THE CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C PWSCYL WILL BE CALLED AGAIN WITH INTL = 1 . C C C DIMENSION F(IDIMF,1) DIMENSION BDA(1) ,BDB(1) ,BDC(1) ,BDD(1) , 1 W(1) C C CHECK FOR INVALID PARAMETERS. C IERROR = 0 IF (A .LT. 0.) IERROR = 1 IF (A .GE. B) IERROR = 2 IF (MBDCND.LE.0 .OR. MBDCND.GE.7) IERROR = 3 IF (C .GE. D) IERROR = 4 IF (NCHECK(N).GT.1 .OR. N.LE.2) IERROR = 5 IF (NBDCND.LE.-1 .OR. NBDCND.GE.5) IERROR = 6 IF (A.EQ.0. .AND. (MBDCND.EQ.3 .OR. MBDCND.EQ.4)) IERROR = 7 IF (A.GT.0. .AND. MBDCND.GE.5) IERROR = 8 IF (A.EQ.0. .AND. ELMBDA.NE.0. .AND. MBDCND.GE.5) IERROR = 9 IF (IDIMF .LT. M+1) IERROR = 10 IF (IERROR .NE. 0) RETURN MP1 = M+1 DELTAR = (B-A)/FLOAT(M) DLRBY2 = DELTAR/2. DLRSQ = DELTAR**2 NP1 = N+1 DELTHT = (D-C)/FLOAT(N) DLTHSQ = DELTHT**2 NP = NBDCND+1 C C DEFINE RANGE OF INDICES I AND J FOR UNKNOWNS U(I,J). C MSTART = 2 MSTOP = M GO TO (104,103,102,101,101,102),MBDCND 101 MSTART = 1 GO TO 104 102 MSTART = 1 103 MSTOP = MP1 104 MUNK = MSTOP-MSTART+1 NSTART = 1 NSTOP = N GO TO (108,105,106,107,108),NP 105 NSTART = 2 GO TO 108 106 NSTART = 2 107 NSTOP = NP1 108 NUNK = NSTOP-NSTART+1 C C DEFINE A,B,C COEFFICIENTS IN W-ARRAY. C ID1 = 6*NUNK+5*MUNK ID2 = ID1+MUNK ID3 = ID2+MUNK ID4 = ID3+MUNK ID5 = ID1-5*MUNK ID6 = ID5+MUNK ISTART = 1 A1 = 2./DLRSQ IJ = 0 IF (MBDCND.EQ.3 .OR. MBDCND.EQ.4) IJ = 1 IF (MBDCND .LE. 4) GO TO 109 W(ID1+1) = 0. W(ID2+1) = -2.*A1 W(ID3+1) = 2.*A1 ISTART = 2 IJ = 1 109 DO 110 I=ISTART,MUNK R = A+FLOAT(I-IJ)*DELTAR J = ID5+I W(J) = R J = ID6+I W(J) = 1./R**2 J = ID1+I W(J) = (R-DLRBY2)/(R*DLRSQ) J = ID3+I W(J) = (R+DLRBY2)/(R*DLRSQ) K = ID6+I J = ID2+I W(J) = -A1+ELMBDA*W(K) 110 CONTINUE GO TO (114,111,112,113,114,112),MBDCND 111 W(ID2) = A1 GO TO 114 112 W(ID2) = A1 113 W(ID3+1) = A1*FLOAT(ISTART) 114 CONTINUE C C ENTER BOUNDARY DATA FOR R-BOUNDARIES. C GO TO (115,115,117,117,119,119),MBDCND 115 A1 = W(ID1+1) DO 116 J=NSTART,NSTOP F(2,J) = F(2,J)-A1*F(1,J) 116 CONTINUE GO TO 119 117 A1 = 2.*DELTAR*W(ID1+1) DO 118 J=NSTART,NSTOP F(1,J) = F(1,J)+A1*BDA(J) 118 CONTINUE 119 GO TO (120,122,122,120,120,122),MBDCND 120 A1 = W(ID4) DO 121 J=NSTART,NSTOP F(M,J) = F(M,J)-A1*F(MP1,J) 121 CONTINUE GO TO 124 122 A1 = 2.*DELTAR*W(ID4) DO 123 J=NSTART,NSTOP F(MP1,J) = F(MP1,J)-A1*BDB(J) 123 CONTINUE C C ENTER BOUNDARY DATA FOR Z-BOUNDARIES. C 124 A1 = 1./DLTHSQ L = ID5-MSTART+1 GO TO (134,125,125,127,127),NP 125 DO 126 I=MSTART,MSTOP F(I,2) = F(I,2)-A1*F(I,1) 126 CONTINUE GO TO 129 127 A1 = 2./DELTHT DO 128 I=MSTART,MSTOP F(I,1) = F(I,1)+A1*BDC(I) 128 CONTINUE 129 A1 = 1./DLTHSQ GO TO (134,130,132,132,130),NP 130 DO 131 I=MSTART,MSTOP F(I,N) = F(I,N)-A1*F(I,NP1) 131 CONTINUE GO TO 134 132 A1 = 2./DELTHT DO 133 I=MSTART,MSTOP F(I,NP1) = F(I,NP1)-A1*BDD(I) 133 CONTINUE 134 CONTINUE C C ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A C SOLUTION. C PERTRB = 0. IF (ELMBDA) 146,136,135 135 IERROR = 11 GO TO 146 136 W(ID5+1) = .5*(W(ID5+2)-DLRBY2) GO TO (146,146,138,146,146,137),MBDCND 137 W(ID5+1) = .5*W(ID5+1) 138 GO TO (140,146,146,139,146),NP 139 A2 = 2. GO TO 141 140 A2 = 1. 141 K = ID5+MUNK W(K) = .5*(W(K-1)+DLRBY2) S = 0. DO 143 I=MSTART,MSTOP S1 = 0. NSP1 = NSTART+1 NSTM1 = NSTOP-1 DO 142 J=NSP1,NSTM1 S1 = S1+F(I,J) 142 CONTINUE K = I+L S = S+(A2*S1+F(I,NSTART)+F(I,NSTOP))*W(K) 143 CONTINUE S2 = FLOAT(M)*A+(.75+FLOAT((M-1)*(M+1)))*DLRBY2 IF (MBDCND .EQ. 3) S2 = S2+.25*DLRBY2 S1 = (2.+A2*FLOAT(NUNK-2))*S2 PERTRB = S/S1 DO 145 I=MSTART,MSTOP DO 144 J=NSTART,NSTOP F(I,J) = F(I,J)-PERTRB 144 CONTINUE 145 CONTINUE 146 CONTINUE C C MULTIPLY I-TH EQUATION THROUGH BY DELTHT**2 TO PUT EQUATION INTO C CORRECT FORM FOR SUBROUTINE POIS. C DO 148 I=MSTART,MSTOP K = I-MSTART+1 J = ID1+K W(J) = W(J)*DLTHSQ J = ID2+K W(J) = W(J)*DLTHSQ J = ID3+K W(J) = W(J)*DLTHSQ DO 147 J=NSTART,NSTOP F(I,J) = F(I,J)*DLTHSQ 147 CONTINUE 148 CONTINUE NPEROD = NBDCND IF (NPEROD .NE. 4) GO TO 152 C C REVERSE COLUMNS SO THE DERIVATIVE SPECIFIED-SOLUTION SPECIFIED C BOUNDARY CONDITIONS BECOME SOLUTION SPECIFIED-DERIVATIVE SPECIFIED C NPEROD = 2 ID5 = 1 K = N/2 149 DO 151 J=1,K L = N+1-J DO 150 I=MSTART,MSTOP S = F(I,L) F(I,L) = F(I,J) F(I,J) = S 150 CONTINUE 151 CONTINUE GO TO (152,156),ID5 C C CALL POIS TO SOLVE THE SYSTEM OF EQUATIONS. C 152 CALL POIS (INTL,NPEROD,NUNK,MBDCND,MUNK,W(ID1+1),W(ID2+1), 1 W(ID3+1),IDIMF,F(MSTART,NSTART),IERR1,W) GO TO (154,156,156,156,153),NP 153 ID5 = 2 GO TO 149 154 DO 155 I=MSTART,MSTOP F(I,NP1) = F(I,1) 155 CONTINUE 156 CONTINUE RETURN END SUBROUTINE PWSSSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND, POI01387 1 BDPS,BDPF,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE PWSSSP SOLVES A FINITE DIFFERENCE APPROXIMATION TO THE C HELMHOLTZ EQUATION IN SPHERICAL COORDINATES AND ON THE SURFACE OF C THE UNIT SPHERE (RADIUS OF 1): C C (1/SIN(THETA))(D/DTHETA)(SIN(THETA)(D/DTHETA)U) C C + (1/SIN(THETA)**2)(D/DPHI)(D/DPHI)U C C + LAMBDA*U = F(THETA,PHI) C C WHERE THETA IS COLATITUDE AND PHI IS LONGITUDE. C C THE ARGUMENTS ARE DEFINED AS: C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C INTL C = 0 ON INITIAL ENTRY TO PWSSSP OR IF N,NBDCND,PS OR PF ARE C CHANGED FROM A PREVIOUS CALL C = 1 IF PS,PF,N AND NBDCND ARE UNCHANGED FROM A PREVIOUS CALL C C NOTE: A CALL WITH INTL = 1 IS ABOUT 1 PERCENT FASTER THAN A C CALL WITH INTL = 0 . C C TS,TF C THE RANGE OF THETA (COLATITUDE), I.E., TS .LE. THETA .LE. TF. C TS MUST BE LESS THAN TF. TS AND TF ARE IN RADIANS. A TS OF C ZERO CORRESPONDS TO THE NORTH POLE AND A TF OF PI CORRESPONDS TO C THE SOUTH POLE. C C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (TS,TF) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C THETA-DIRECTION GIVEN BY THETA(I) = (I-1)DTHETA+TS FOR C I = 1,2,...,M+1, WHERE DTHETA = (TF-TS)/M IS THE PANEL WIDTH. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION AT THETA = TS AND C THETA = TF. C C = 1 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THETA = TF. C = 2 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW). C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS AND THETA = TF (SEE NOTES 1,2 C BELOW). C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE C SOLUTION IS SPECIFIED AT THETA = TF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE C SOLUTION IS SPECIFIED AT THETA = TF. C = 6 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW). C = 7 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE C SOLUTION IS UNSPECIFIED AT THETA = TF = PI. C = 8 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE C SOLUTION IS UNSPECIFIED AT THETA = TF = PI. C = 9 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND C THETA = TF = PI. C C NOTES: 1. IF TS = 0, DO NOT USE MBDCND = 3,4, OR 8, BUT C INSTEAD USE MBDCND = 5,6, OR 9 . C 2. IF TF = PI, DO NOT USE MBDCND = 2,3, OR 6, BUT C INSTEAD USE MBDCND = 7,8, OR 9 . C C BDTS C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = TS. WHEN MBDCND = 3,4, OR 8, C C BDTS(J) = (D/DTHETA)U(TS,PHI(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTS IS A DUMMY VARIABLE. C C BDTF C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = TF. WHEN MBDCND = 2,3, OR 6, C C BDTF(J) = (D/DTHETA)U(TF,PHI(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTF IS A DUMMY VARIABLE. C C PS,PF C THE RANGE OF PHI (LONGITUDE), I.E., PS .LE. PHI .LE. PF. PS C MUST BE LESS THAN PF. PS AND PF ARE IN RADIANS. IF PS = 0 AND C PF = 2*PI, PERIODIC BOUNDARY CONDITIONS ARE USUALLY PRESCRIBED. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (PS,PF) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C PHI-DIRECTION GIVEN BY PHI(J) = (J-1)DPHI+PS FOR C J = 1,2,...,N+1, WHERE DPHI = (PF-PS)/N IS THE PANEL WIDTH. N C MUST BE OF THE FORM (2**P)(3**Q)(5**R) WHERE P, Q, AND R ARE ANY C NON-NEGATIVE INTEGERS. N MUST BE GREATER THAN 2 . C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION AT PHI = PS AND C PHI = PF. C C = 0 IF THE SOLUTION IS PERIODIC IN PHI, I.E., C U(I,1) = U(I,N+1). C = 1 IF THE SOLUTION IS SPECIFIED AT PHI = PS AND PHI = PF C (SEE NOTE BELOW). C = 2 IF THE SOLUTION IS SPECIFIED AT PHI = PS (SEE NOTE BELOW) C AND THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS C SPECIFIED AT PHI = PF. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS C SPECIFIED AT PHI = PS AND PHI = PF. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI IS C SPECIFIED AT PS AND THE SOLUTION IS SPECIFIED AT PHI = PF C (SEE NOTE BELOW). C C NOTE: NBDCND = 1,2, OR 4 CANNOT BE USED WITH C MBDCND = 5,6,7,8, OR 9 (THE FORMER INDICATES THAT THE C SOLUTION IS SPECIFIED AT A POLE, THE LATTER C INDICATES THAT THE SOLUTION IS UNSPECIFIED). C USE INSTEAD C MBDCND = 1 OR 2 . C C BDPS C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI AT C PHI = PS. WHEN NBDCND = 3 OR 4, C C BDPS(I) = (D/DPHI)U(THETA(I),PS), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDPS IS A DUMMY VARIABLE. C C BDPF C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO PHI AT C PHI = PF. WHEN NBDCND = 2 OR 3, C C BDPF(I) = (D/DPHI)U(THETA(I),PF), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDPF IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, PWSSSP WILL C ATTEMPT TO FIND A SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUE OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY). C FOR I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(THETA(I),PHI(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ ------------ ------------ C C 1 U(TS,PHI(J)) U(TF,PHI(J)) C 2 U(TS,PHI(J)) F(TF,PHI(J)) C 3 F(TS,PHI(J)) F(TF,PHI(J)) C 4 F(TS,PHI(J)) U(TF,PHI(J)) C 5 F(0,PS) U(TF,PHI(J)) J = 1,2,...,N+1 C 6 F(0,PS) F(TF,PHI(J)) C 7 U(TS,PHI(J)) F(PI,PS) C 8 F(TS,PHI(J)) F(PI,PS) C 9 F(0,PS) F(PI,PS) C C NBDCND F(I,1) F(I,N+1) C ------ -------------- -------------- C C 0 F(THETA(I),PS) F(THETA(I),PS) C 1 U(THETA(I),PS) U(THETA(I),PF) C 2 U(THETA(I),PS) F(THETA(I),PF) I = 1,2,...,M+1 C 3 F(THETA(I),PS) F(THETA(I),PF) C 4 F(THETA(I),PS) U(THETA(I),PF) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING PWSSSP. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 11(M+1)+6(N+1). C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (THETA(I),PHI(J)), C I = 1,2,...,M+1, J = 1,2,...,N+1 . C C PERTRB C IF ONE SPECIFIES A COMBINATION OF PERIODIC, DERIVATIVE OR C UNSPECIFIED BOUNDARY CONDITIONS FOR A POISSON EQUATION C (LAMBDA = 0), A SOLUTION MAY NOT EXIST. PERTRB IS A CONSTANT, C CALCULATED AND SUBTRACTED FROM F, WHICH ENSURES THAT A SOLUTION C EXISTS. PWSSSP THEN COMPUTES THIS SOLUTION, WHICH IS A LEAST C SQUARES SOLUTION TO THE ORIGINAL APPROXIMATION. THIS SOLUTION C IS NOT UNIQUE AND IS UNNORMALIZED. THE VALUE OF PERTRB SHOULD C BE SMALL COMPARED TO THE RIGHT SIDE F. OTHERWISE , A SOLUTION C IS OBTAINED TO AN ESSENTIALLY DIFFERENT PROBLEM. THIS COMPARISON C SHOULD ALWAYS BE MADE TO INSURE THAT A MEANINGFUL SOLUTION HAS C BEEN OBTAINED C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 8, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 TS .LT. 0 OR TF .GT. PI C = 2 TS .GE. TF. C = 3 MBDCND .LT. 1 OR MBDCND .GT. 9 . C = 4 PS .GE. PF. C = 5 N IS NOT OF THE FORM (2**P)(3**Q)(5**R) OR N .LE. 2 . C = 6 NBDCND .LT. 0 OR NBDCND .GT. 4 . C = 7 AN NBDCND OF 1,2, OR 4 IS USED WITH AN C MBDCND OF 5,6,7,8, OR 9 . C = 8 ELMBDA .GT. 0 . C = 9 IDIMF .LT. M+1 . C = 10 TS=0 AND MBDCND=3,4 OR 8 OR TF=PI AND MBDCND=2,3 OR 6 C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSIBLY INCORRECT C CALL TO PWSSSP, THE USER SHOULD TEST IERROR AFTER A CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C PWSSSP WILL BE CALLED AGAIN WITH INTL = 1 . C C C DIMENSION F(IDIMF,1) ,BDTS(1) ,BDTF(1) ,BDPS(1) , 1 BDPF(1) ,W(1) NBR = NBDCND+1 PI = 3.14159265358979 EPS = APXEPS(DUM) IERROR = 9 IF (IDIMF-M) 119,119,101 101 IERROR = 1 IF (TS) 119,102,102 102 IF (PI-TF+EPS) 119,119,103 103 IERROR = 2 IF (TF-TS) 119,119,104 104 IERROR = 3 IF (MBDCND) 119,119,105 105 IF (MBDCND-9) 106,106,119 106 IERROR = 4 IF (PF-PS) 119,119,107 107 IERROR = 5 IF (NCHECK(N)-2) 108,119,108 108 IERROR = 6 IF (NBDCND) 119,109,109 109 IF (NBDCND-4) 110,110,119 110 IERROR = 7 IF (MBDCND-4) 112,112,111 111 GO TO (112,119,119,112,119),NBR 112 IERROR = 8 IF (ELMBDA) 113,113,118 113 IERROR = 10 IF (TS) 114,114,115 114 GO TO (115,115,118,118,115,115,115,118,115),MBDCND 115 IF (PI-TF-EPS) 116,116,117 116 GO TO (117,118,118,117,117,118,117,117,117),MBDCND 117 IERROR = 0 118 MP1 = M+1 IW1 = 6*(N+1)+5*MP1+1 IW2 = IW1+MP1 IW3 = IW2+MP1 IW4 = IW3+MP1 IW5 = IW4+MP1 IW6 = IW5+MP1 CALL PWSSS1 (INTL,TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS, 1 BDPF,ELMBDA,F,IDIMF,PERTRB,W(IW1),W(IW2),W(IW3), 2 W(IW4),W(IW5),W(IW6),W) 119 RETURN END SUBROUTINE PWSSS1 (INTL,TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND, POI01701 1 BDPS,BDPF,ELMBDA,F,IDIMF,PERTRB,AM,BM,CM,SN, 2 SS,SINT,D) DIMENSION F(IDIMF,1) ,BDTS(1) ,BDTF(1) ,BDPS(1) , 1 BDPF(1) ,AM(1) ,BM(1) ,CM(1) , 2 SS(1) ,SN(1) ,D(1) ,SINT(1) C C MACHINE DEPENDENT CONSTANT C C PI=3.1415926535897932384626433832795028841971693993751058209749446 C PI = 3.14159265358979 TPI = PI+PI HPI = PI/2. MP1 = M+1 NP1 = N+1 FN = N FM = M DTH = (TF-TS)/FM HDTH = DTH/2. TDT = DTH+DTH DPHI = (PF-PS)/FN TDP = DPHI+DPHI DPHI2 = DPHI*DPHI EDP2 = ELMBDA*DPHI2 DTH2 = DTH*DTH CP = 4./(FN*DTH2) WP = FN*SIN(HDTH)/4. DO 102 I=1,MP1 FIM1 = I-1 THETA = FIM1*DTH+TS SINT(I) = SIN(THETA) IF (SINT(I)) 101,102,101 101 T1 = 1./(DTH2*SINT(I)) AM(I) = T1*SIN(THETA-HDTH) CM(I) = T1*SIN(THETA+HDTH) BM(I) = -AM(I)-CM(I)+ELMBDA 102 CONTINUE INP = 0 ISP = 0 C C BOUNDARY CONDITION AT THETA=TS C MBR = MBDCND+1 GO TO (103,104,104,105,105,106,106,104,105,106),MBR 103 ITS = 1 GO TO 107 104 AT = AM(2) ITS = 2 GO TO 107 105 AT = AM(1) ITS = 1 CM(1) = AM(1)+CM(1) GO TO 107 106 AT = AM(2) INP = 1 ITS = 2 C C BOUNDARY CONDITION THETA=TF C 107 GO TO (108,109,110,110,109,109,110,111,111,111),MBR 108 ITF = M GO TO 112 109 CT = CM(M) ITF = M GO TO 112 110 CT = CM(M+1) AM(M+1) = AM(M+1)+CM(M+1) ITF = M+1 GO TO 112 111 ITF = M ISP = 1 CT = CM(M) C C COMPUTE HOMOGENEOUS SOLUTION WITH SOLUTION AT POLE EQUAL TO ONE C 112 ITSP = ITS+1 ITFM = ITF-1 WTS = SINT(ITS+1)*AM(ITS+1)/CM(ITS) WTF = SINT(ITF-1)*CM(ITF-1)/AM(ITF) MUNK = ITF-ITS+1 IF (ISP) 116,116,113 113 D(ITS) = CM(ITS)/BM(ITS) DO 114 I=ITSP,M D(I) = CM(I)/(BM(I)-AM(I)*D(I-1)) 114 CONTINUE SS(M) = -D(M) IID = M-ITS DO 115 II=1,IID I = M-II SS(I) = -D(I)*SS(I+1) 115 CONTINUE SS(M+1) = 1. 116 IF (INP) 120,120,117 117 SN(1) = 1. D(ITF) = AM(ITF)/BM(ITF) IID = ITF-2 DO 118 II=1,IID I = ITF-II D(I) = AM(I)/(BM(I)-CM(I)*D(I+1)) 118 CONTINUE SN(2) = -D(2) DO 119 I=3,ITF SN(I) = -D(I)*SN(I-1) 119 CONTINUE C C BOUNDARY CONDITIONS AT PHI=PS C 120 NBR = NBDCND+1 WPS = 1. WPF = 1. GO TO (121,122,122,123,123),NBR 121 JPS = 1 GO TO 124 122 JPS = 2 GO TO 124 123 JPS = 1 WPS = .5 C C BOUNDARY CONDITION AT PHI=PF C 124 GO TO (125,126,127,127,126),NBR 125 JPF = N GO TO 128 126 JPF = N GO TO 128 127 WPF = .5 JPF = N+1 128 JPSP = JPS+1 JPFM = JPF-1 NUNK = JPF-JPS+1 FJJ = JPFM-JPSP+1 C C SCALE COEFFICIENTS FOR SUBROUTINE POIS C DO 129 I=ITS,ITF CF = DPHI2*SINT(I)*SINT(I) AM(I) = CF*AM(I) BM(I) = CF*BM(I) CM(I) = CF*CM(I) 129 CONTINUE ISING = 0 GO TO (130,138,138,130,138,138,130,138,130,130),MBR 130 GO TO (131,138,138,131,138),NBR 131 IF (ELMBDA) 138,132,132 132 ISING = 1 SUM = WTS*WPS+WTS*WPF+WTF*WPS+WTF*WPF IF (INP) 134,134,133 133 SUM = SUM+WP 134 IF (ISP) 136,136,135 135 SUM = SUM+WP 136 SUM1 = 0. DO 137 I=ITSP,ITFM SUM1 = SUM1+SINT(I) 137 CONTINUE SUM = SUM+FJJ*(SUM1+WTS+WTF) SUM = SUM+(WPS+WPF)*SUM1 HNE = SUM 138 GO TO (146,142,142,144,144,139,139,142,144,139),MBR 139 IF (NBDCND-3) 146,140,146 140 YHLD = F(1,JPS)-4./(FN*DPHI*DTH2)*(BDPF(2)-BDPS(2)) DO 141 J=1,NP1 F(1,J) = YHLD 141 CONTINUE GO TO 146 142 DO 143 J=JPS,JPF F(2,J) = F(2,J)-AT*F(1,J) 143 CONTINUE GO TO 146 144 DO 145 J=JPS,JPF F(1,J) = F(1,J)+TDT*BDTS(J)*AT 145 CONTINUE 146 GO TO (154,150,152,152,150,150,152,147,147,147),MBR 147 IF (NBDCND-3) 154,148,154 148 YHLD = F(M+1,JPS)-4./(FN*DPHI*DTH2)*(BDPF(M)-BDPS(M)) DO 149 J=1,NP1 F(M+1,J) = YHLD 149 CONTINUE GO TO 154 150 DO 151 J=JPS,JPF F(M,J) = F(M,J)-CT*F(M+1,J) 151 CONTINUE GO TO 154 152 DO 153 J=JPS,JPF F(M+1,J) = F(M+1,J)-TDT*BDTF(J)*CT 153 CONTINUE 154 GO TO (159,155,155,157,157),NBR 155 DO 156 I=ITS,ITF F(I,2) = F(I,2)-F(I,1)/(DPHI2*SINT(I)*SINT(I)) 156 CONTINUE GO TO 159 157 DO 158 I=ITS,ITF F(I,1) = F(I,1)+TDP*BDPS(I)/(DPHI2*SINT(I)*SINT(I)) 158 CONTINUE 159 GO TO (164,160,162,162,160),NBR 160 DO 161 I=ITS,ITF F(I,N) = F(I,N)-F(I,N+1)/(DPHI2*SINT(I)*SINT(I)) 161 CONTINUE GO TO 164 162 DO 163 I=ITS,ITF F(I,N+1) = F(I,N+1)-TDP*BDPF(I)/(DPHI2*SINT(I)*SINT(I)) 163 CONTINUE 164 CONTINUE PERTRB = 0. IF (ISING) 165,176,165 165 SUM = WTS*WPS*F(ITS,JPS)+WTS*WPF*F(ITS,JPF)+WTF*WPS*F(ITF,JPS)+ 1 WTF*WPF*F(ITF,JPF) IF (INP) 167,167,166 166 SUM = SUM+WP*F(1,JPS) 167 IF (ISP) 169,169,168 168 SUM = SUM+WP*F(M+1,JPS) 169 DO 171 I=ITSP,ITFM SUM1 = 0. DO 170 J=JPSP,JPFM SUM1 = SUM1+F(I,J) 170 CONTINUE SUM = SUM+SINT(I)*SUM1 171 CONTINUE SUM1 = 0. SUM2 = 0. DO 172 J=JPSP,JPFM SUM1 = SUM1+F(ITS,J) SUM2 = SUM2+F(ITF,J) 172 CONTINUE SUM = SUM+WTS*SUM1+WTF*SUM2 SUM1 = 0. SUM2 = 0. DO 173 I=ITSP,ITFM SUM1 = SUM1+SINT(I)*F(I,JPS) SUM2 = SUM2+SINT(I)*F(I,JPF) 173 CONTINUE SUM = SUM+WPS*SUM1+WPF*SUM2 PERTRB = SUM/HNE DO 175 J=1,NP1 DO 174 I=1,MP1 F(I,J) = F(I,J)-PERTRB 174 CONTINUE 175 CONTINUE C C SCALE RIGHT SIDE FOR SUBROUTINE POIS C 176 DO 178 I=ITS,ITF CF = DPHI2*SINT(I)*SINT(I) DO 177 J=JPS,JPF F(I,J) = CF*F(I,J) 177 CONTINUE 178 CONTINUE NBDC = NBDCND IF (NBDCND-4) 182,179,182 179 JFN = NUNK/2 DO 181 J=1,JFN JFRD = J+JPS-1 JBRD = JPF-J+1 DO 180 I=ITS,ITF HLD = F(I,JFRD) F(I,JFRD) = F(I,JBRD) F(I,JBRD) = HLD 180 CONTINUE 181 CONTINUE NBDC = 2 182 CALL POIS (INTL,NBDC,NUNK,MBDCND,MUNK,AM(ITS),BM(ITS),CM(ITS), 1 IDIMF,F(ITS,JPS),IERROR,D) IF (NBDCND-4) 186,183,186 183 DO 185 J=1,JFN JFRD = J+JPS-1 JBRD = JPF-J+1 DO 184 I=ITS,ITF HLD = F(I,JFRD) F(I,JFRD) = F(I,JBRD) F(I,JBRD) = HLD 184 CONTINUE 185 CONTINUE 186 IF (ISING) 194,194,187 187 IF (INP) 191,191,188 188 IF (ISP) 189,189,194 189 DO 190 J=1,NP1 F(1,J) = 0. 190 CONTINUE GO TO 217 191 IF (ISP) 194,194,192 192 DO 193 J=1,NP1 F(M+1,J) = 0. 193 CONTINUE GO TO 217 194 IF (INP) 201,201,195 195 SUM = WPS*F(ITS,JPS)+WPF*F(ITS,JPF) DO 196 J=JPSP,JPFM SUM = SUM+F(ITS,J) 196 CONTINUE DFN = CP*SUM DNN = CP*((WPS+WPF+FJJ)*(SN(2)-1.))+ELMBDA DSN = CP*(WPS+WPF+FJJ)*SN(M) IF (ISP) 197,197,202 197 CNP = (F(1,1)-DFN)/DNN DO 199 I=ITS,ITF HLD = CNP*SN(I) DO 198 J=JPS,JPF F(I,J) = F(I,J)+HLD 198 CONTINUE 199 CONTINUE DO 200 J=1,NP1 F(1,J) = CNP 200 CONTINUE GO TO 217 201 IF (ISP) 217,217,202 202 SUM = WPS*F(ITF,JPS)+WPF*F(ITF,JPF) DO 203 J=JPSP,JPFM SUM = SUM+F(ITF,J) 203 CONTINUE DFS = CP*SUM DSS = CP*((WPS+WPF+FJJ)*(SS(M)-1.))+ELMBDA DNS = CP*(WPS+WPF+FJJ)*SS(2) IF (INP) 204,204,208 204 CSP = (F(M+1,1)-DFS)/DSS DO 206 I=ITS,ITF HLD = CSP*SS(I) DO 205 J=JPS,JPF F(I,J) = F(I,J)+HLD 205 CONTINUE 206 CONTINUE DO 207 J=1,NP1 F(M+1,J) = CSP 207 CONTINUE GO TO 217 208 RTN = F(1,1)-DFN RTS = F(M+1,1)-DFS IF (ISING) 210,210,209 209 CSP = 0. CNP = RTN/DNN GO TO 213 210 IF (ABS(DNN)-ABS(DSN)) 212,212,211 211 DEN = DSS-DNS*DSN/DNN RTS = RTS-RTN*DSN/DNN CSP = RTS/DEN CNP = (RTN-CSP*DNS)/DNN GO TO 213 212 DEN = DNS-DSS*DNN/DSN RTN = RTN-RTS*DNN/DSN CSP = RTN/DEN CNP = (RTS-DSS*CSP)/DSN 213 DO 215 I=ITS,ITF HLD = CNP*SN(I)+CSP*SS(I) DO 214 J=JPS,JPF F(I,J) = F(I,J)+HLD 214 CONTINUE 215 CONTINUE DO 216 J=1,NP1 F(1,J) = CNP F(M+1,J) = CSP 216 CONTINUE 217 IF (NBDCND) 220,218,220 218 DO 219 I=1,MP1 F(I,JPF+1) = F(I,JPS) 219 CONTINUE 220 RETURN END SUBROUTINE PWSCSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND, POI02058 1 BDRS,BDRF,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE PWSCSP SOLVES A FINITE DIFFERENCE APPROXIMATION TO THE C MODIFIED HELMHOLTZ EQUATION IN SPHERICAL COORDINATES ASSUMING C AXISYMMETRY (NO DEPENDENCE ON LONGITUDE) C C (1/R**2)(D/DR)((R**2)(D/DR)U) C C + (1/(R**2)SIN(THETA))(D/DTHETA)(SIN(THETA)(D/DTHETA)U) C C + (LAMBDA/(RSIN(THETA))**2)U = F(THETA,R). C THIS TWO DIMENSIONAL MODIFIED HELMHOLTZ EQUATION RESULTS FROM C THE FOURIER TRANSFORM OF THE THREE DIMENSIONAL POISSON EQUATION C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C INTL C = 0 ON INITIAL ENTRY TO PWSCSP OR IF ANY OF THE ARGUMENTS C RS, RF, N, NBDCND ARE CHANGED FROM A PREVIOUS CALL. C = 1 IF RS, RF, N, NBDCND ARE ALL UNCHANGED FROM PREVIOUS CALL C TO PWSCSP. C C NOTE A CALL WITH INTL=0 TAKES APPROXIMATELY 1.5 TIMES AS C MUCH TIME AS A CALL WITH INTL = 1 . ONCE A CALL WITH C INTL = 0 HAS BEEN MADE THEN SUBSEQUENT SOLUTIONS C CORRESPONDING TO DIFFERENT F, BDTS, BDTF, BDRS, BDRF CAN C BE OBTAINED FASTER WITH INTL = 1 SINCE INITIALIZATION IS C NOT REPEATED. C C TS,TF C THE RANGE OF THETA (COLATITUDE), I.E., TS .LE. THETA .LE. TF. C TS MUST BE LESS THAN TF. TS AND TF ARE IN RADIANS. A TS OF C ZERO CORRESPONDS TO THE NORTH POLE AND A TF OF PI CORRESPONDS C TO THE SOUTH POLE. C M C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (TS,TF) IS C SUBDIVIDED. HENCE, THERE WILL BE M+1 GRID POINTS IN THE C THETA-DIRECTION GIVEN BY THETA(K) = (I-1)DTHETA+TS FOR C I = 1,2,...,M+1, WHERE DTHETA = (TF-TS)/M IS THE PANEL WIDTH. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION AT THETA = TS AND C THETA = TF. C C = 1 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THETA = TF. C = 2 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW). C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS AND THETA = TF (SEE NOTES 1,2 C BELOW). C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE C SOLUTION IS SPECIFIED AT THETA = TF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE C SOLUTION IS SPECIFIED AT THETA = TF. C = 6 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND THE C DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TF (SEE NOTE 2 BELOW). C = 7 IF THE SOLUTION IS SPECIFIED AT THETA = TS AND THE C SOLUTION IS UNSPECIFIED AT THETA = TF = PI. C = 8 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TS (SEE NOTE 1 BELOW) AND THE SOLUTION C IS UNSPECIFIED AT THETA = TF = PI. C = 9 IF THE SOLUTION IS UNSPECIFIED AT THETA = TS = 0 AND C THETA = TF = PI. C C NOTES: 1. IF TS = 0, DO NOT USE MBDCND = 3,4, OR 8, BUT C INSTEAD USE MBDCND = 5,6, OR 9 . C 2. IF TF = PI, DO NOT USE MBDCND = 2,3, OR 6, BUT C INSTEAD USE MBDCND = 7,8, OR 9 . C C BDTS C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = TS. WHEN MBDCND = 3,4, OR 8, C C BDTS(J) = (D/DTHETA)U(TS,R(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTS IS A DUMMY VARIABLE. C C BDTF C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO THETA AT C THETA = TF. WHEN MBDCND = 2,3, OR 6, C C BDTF(J) = (D/DTHETA)U(TF,R(J)), J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTF IS A DUMMY VARIABLE. C C RS,RF C THE RANGE OF R, I.E., RS .LE. R .LT. RF. RS MUST BE LESS THAN C RF. RS MUST BE NON-NEGATIVE. C C N C THE NUMBER OF PANELS INTO WHICH THE INTERVAL (RS,RF) IS C SUBDIVIDED. HENCE, THERE WILL BE N+1 GRID POINTS IN THE C R-DIRECTION GIVEN BY R(J) = (J-1)DR+RS FOR J = 1,2,...,N+1, C WHERE DR = (RF-RS)/N IS THE PANEL WIDTH. LET K BE AN INTEGER C GREATER THAN ONE, THEN N MUST HAVE THE FOLLOWING FORM DEPENDING C ON NBDCND (SEE BELOW). C C IF NBDCND = 2,4, OR 6, N MUST HAVE THE FORM 2**K-1 . C IF NBDCND = 1 OR 5, N MUST HAVE THE FORM 2**K. C IF NBDCND = 3, N MUST HAVE THE FORM 2**K-2 . C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION AT R = RS AND R = RF. C C = 1 IF THE SOLUTION IS SPECIFIED AT R = RS AND R = RF. C = 2 IF THE SOLUTION IS SPECIFIED AT R = RS AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO R IS SPECIFIED AT R = RF. C = 3 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = RS AND R = RF. C = 4 IF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R IS C SPECIFIED AT RS AND THE SOLUTION IS SPECIFIED AT R = RF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT R = RS = 0 (SEE NOTE C BELOW) AND THE SOLUTION IS SPECIFIED AT R = RF. C = 6 IF THE SOLUTION IS UNSPECIFIED AT R = RS = 0 (SEE NOTE C BELOW) AND THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO C R IS SPECIFIED AT R = RF. C C NOTE: NBDCND = 5 OR 6 CANNOT BE USED WITH C MBDCND = 1,2,4,5, OR 7 (THE FORMER INDICATES THAT THE C SOLUTION IS UNSPECIFIED AT R = 0, THE LATTER C INDICATES THAT THE SOLUTION IS SPECIFIED). C USE INSTEAD C NBDCND = 1 OR 2 . C C BDRS C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = RS. C WHEN NBDCND = 3 OR 4, C C BDRS(I) = (D/DR)U(THETA(I),RS), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDRS IS A DUMMY VARIABLE. C C BDRF C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT SPECIFIES THE VALUES C OF THE DERIVATIVE OF THE SOLUTION WITH RESPECT TO R AT R = RF. C WHEN NBDCND = 2,3, OR 6, C C BDRF(I) = (D/DR)U(THETA(I),RF), I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDRF IS A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ EQUATION. IF C LAMBDA .GT. 0, A SOLUTION MAY NOT EXIST. HOWEVER, PWSCSP WILL C ATTEMPT TO FIND A SOLUTION. IF NBDCND = 5 OR 6 OR C MBDCND = 5,6,7,8, OR 9, ELMBDA MUST BE ZERO. C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUE OF THE RIGHT C SIDE OF THE HELMHOLTZ EQUATION AND BOUNDARY VALUES (IF ANY). C FOR I = 2,3,...,M AND J = 2,3,...,N C C F(I,J) = F(THETA(I),R(J)). C C ON THE BOUNDARIES F IS DEFINED BY C C MBDCND F(1,J) F(M+1,J) C ------ ---------- ---------- C C 1 U(TS,R(J)) U(TF,R(J)) C 2 U(TS,R(J)) F(TF,R(J)) C 3 F(TS,R(J)) F(TF,R(J)) C 4 F(TS,R(J)) U(TF,R(J)) C 5 F(0,R(J)) U(TF,R(J)) J = 1,2,...,N+1 C 6 F(0,R(J)) F(TF,R(J)) C 7 U(TS,R(J)) F(PI,R(J)) C 8 F(TS,R(J)) F(PI,R(J)) C 9 F(0,R(J)) F(PI,R(J)) C C NBDCND F(I,1) F(I,N+1) C ------ -------------- -------------- C C 1 U(THETA(I),RS) U(THETA(I),RF) C 2 U(THETA(I),RS) F(THETA(I),RF) C 3 F(THETA(I),RS) F(THETA(I),RF) C 4 F(THETA(I),RS) U(THETA(I),RF) I = 1,2,...,M+1 C 5 F(TS,0) U(THETA(I),RF) C 6 F(TS,0) F(THETA(I),RF) C C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C NOTE C C IF THE TABLE CALLS FOR BOTH THE SOLUTION U AND THE RIGHT SIDE F C AT A CORNER THEN THE SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY F AS IT APPEARS IN THE C PROGRAM CALLING PWSCSP. THIS PARAMETER IS USED TO SPECIFY THE C VARIABLE DIMENSION OF F. IDIMF MUST BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. LET L = 2**K (SEE THE DEFINITION OF N). THEN THE C LENGTH OF W MUST BE AT LEAST C C 2(L+1)(K-1)+6(M+N)+MAX(4N,6M)+14 C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE DIFFERENCE C APPROXIMATION FOR THE GRID POINT (THETA(I),R(J)), C I = 1,2,...,M+1, J = 1,2,...,N+1 . C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE BOUNDARY CONDITIONS C IS SPECIFIED FOR A POISSON EQUATION (LAMBDA = 0), A SOLUTION MAY C NOT EXIST. PERTRB IS A CONSTANT, CALCULATED AND SUBTRACTED FROM C F, WHICH ENSURES THAT A SOLUTION EXISTS. PWSCSP THEN COMPUTES C THIS SOLUTION, WHICH IS A LEAST SQUARES SOLUTION TO THE ORIGINAL C APPROXIMATION. THIS SOLUTION IS NOT UNIQUE AND IS UNNORMALIZED. C THE VALUE OF PERTRB SHOULD BE SMALL COMPARED TO THE RIGHT SIDE C F. OTHERWISE , A SOLUTION IS OBTAINED TO AN ESSENTIALLY C DIFFERENT PROBLEM. THIS COMPARISON SHOULD ALWAYS BE MADE TO C INSURE THAT A MEANINGFUL SOLUTION HAS BEEN OBTAINED. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBERS 0 AND 10, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 TS .LT. 0 OR TF .GT. PI C = 2 TS .GE. TF. C = 3 MBDCND .LT. 1 OR MBDCND .GT. 9 . C = 4 RS .LT. 0 . C = 5 RS .GE. RF. C = 6 NBDCND .LT. 1 OR NBDCND .GT. 6 . C = 7 N IS NOT OF THE PROPER FORM. C = 8 AN NBDCND OF 5 OR 6 IS USED WITH AN C MBDCND OF 1,2,4,5, OR 7 . C = 9 ELMBDA IS NON-ZERO AND EITHER C MBDCND = 5,6,7,8, OR 9 OR NBDCND = 5 OR 6 . C = 10 ELMBDA .GT. 0 . C = 11 IDIMF .LT. M+1 . C = 12 TS=0 AND MBDCND=3,4 OR 8 OR TF=PI AND MBDCND=2,3 OR 6 C C SINCE THIS IS THE ONLY MEANS OF INDICATING A POSSLIBY INCORRECT C CALL TO PWSCSP, THE USER SHOULD TEST IERROR AFTER A CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C PWSCSP WILL BE CALLED AGAIN WITH INTL = 1 . C C C DIMENSION F(IDIMF,1) ,BDTS(1) ,BDTF(1) ,BDRS(1) , 1 BDRF(1) ,W(1) PI = 3.14159265358979 EPS = APXEPS(DUM) IERROR = 11 IF (IDIMF-M) 132,132,101 101 IERROR = 1 IF (TS) 132,102,102 102 IF (PI-TF+EPS) 132,132,103 103 IERROR = 2 IF (TF-TS) 132,132,104 104 IERROR = 3 IF (MBDCND) 132,132,105 105 IF (MBDCND-9) 106,106,132 106 IERROR = 4 IF (RS) 132,107,107 107 IERROR = 5 IF (RF-RS) 132,132,108 108 IERROR = 6 IF (NBDCND) 132,132,109 109 IF (NBDCND-6) 110,110,132 110 IERROR = 7 NCK = N GO TO (113,111,112,111,113,111),NBDCND 111 NCK = N+1 GO TO 113 112 NCK = N+2 113 K = -1 114 IF (NCK) 116,116,115 115 NCK = NCK/2 K = K+1 GO TO 114 116 GO TO (117,118,119,118,117,118),NBDCND 117 IF (2**K-N) 132,120,132 118 IF (2**K-1-N) 132,120,132 119 IF (2**K-2-N) 132,120,132 120 IERROR = 8 GO TO (121,121,122,121,121,122,121,122,122),MBDCND 121 IF (NBDCND-4) 122,122,132 122 IERROR = 9 IF (ELMBDA) 123,125,123 123 IF (MBDCND-4) 124,124,132 124 IF (NBDCND-4) 125,125,132 125 IERROR = 10 IF (ELMBDA) 126,126,131 126 IERROR = 12 IF (TS) 127,127,130 127 GO TO (128,128,131,131,128,128,128,131,128),MBDCND 128 IF (PI-TF-EPS) 129,129,130 129 GO TO (130,131,131,130,130,131,130,130,130),MBDCND 130 IERROR = 0 131 I1 = 6*M+N+8 I2 = I1+N+1 I3 = I2+N+1 I4 = I3+N+1 CALL PWSCS1 (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND,BDRS, 1 BDRF,ELMBDA,F,IDIMF,PERTRB,W,W(M+2),W(2*M+3), 2 W(3*M+4),W(4*M+5),W(5*M+6),W(6*M+7),W(I1),W(I2), 3 W(I3),W(I4)) 132 RETURN END SUBROUTINE PWSCS1 (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND, POI02399 1 BDRS,BDRF,ELMBDA,F,IDIMF,PERTRB,AM,BM,CM,S, 2 SINT,BMH,AN,BN,CN,R,W) DIMENSION F(IDIMF,1) ,BDRS(1) ,BDRF(1) ,BDTS(1) , 1 BDTF(1) ,AM(1) ,BM(1) ,CM(1) , 2 AN(1) ,BN(1) ,CN(1) ,S(1) , 3 R(1) ,SINT(1) ,BMH(1) ,W(1) MP1 = M+1 DTH = (TF-TS)/FLOAT(M) TDT = DTH+DTH HDTH = DTH/2. SDTS = 1./(DTH*DTH) DO 102 I=1,MP1 THETA = TS+FLOAT(I-1)*DTH SINT(I) = SIN(THETA) IF (SINT(I)) 101,102,101 101 T1 = SDTS/SINT(I) AM(I) = T1*SIN(THETA-HDTH) CM(I) = T1*SIN(THETA+HDTH) BM(I) = -(AM(I)+CM(I)) 102 CONTINUE NP1 = N+1 DR = (RF-RS)/FLOAT(N) HDR = DR/2. TDR = DR+DR DR2 = DR*DR CZR = 6.*DTH/(DR2*(COS(TS)-COS(TF))) DO 103 J=1,NP1 R(J) = RS+FLOAT(J-1)*DR AN(J) = (R(J)-HDR)**2/DR2 CN(J) = (R(J)+HDR)**2/DR2 BN(J) = -(AN(J)+CN(J)) 103 CONTINUE MP = 1 NP = 1 C C BOUNDARY CONDITION AT PHI=PS C GO TO (104,104,105,105,106,106,104,105,106),MBDCND 104 AT = AM(2) ITS = 2 GO TO 107 105 AT = AM(1) ITS = 1 CM(1) = CM(1)+AM(1) GO TO 107 106 ITS = 1 BM(1) = -4.*SDTS CM(1) = -BM(1) C C BOUNDARY CONDITION AT PHI=PF C 107 GO TO (108,109,109,108,108,109,110,110,110),MBDCND 108 CT = CM(M) ITF = M GO TO 111 109 CT = CM(M+1) AM(M+1) = AM(M+1)+CM(M+1) ITF = M+1 GO TO 111 110 ITF = M+1 AM(M+1) = 4.*SDTS BM(M+1) = -AM(M+1) 111 WTS = SINT(ITS+1)*AM(ITS+1)/CM(ITS) WTF = SINT(ITF-1)*CM(ITF-1)/AM(ITF) ITSP = ITS+1 ITFM = ITF-1 C C BOUNDARY CONDITION AT R=RS C ICTR = 0 GO TO (112,112,113,113,114,114),NBDCND 112 AR = AN(2) JRS = 2 GO TO 118 113 AR = AN(1) JRS = 1 CN(1) = CN(1)+AN(1) GO TO 118 114 JRS = 2 ICTR = 1 S(N) = AN(N)/BN(N) DO 115 J=3,N L = N-J+2 S(L) = AN(L)/(BN(L)-CN(L)*S(L+1)) 115 CONTINUE S(2) = -S(2) DO 116 J=3,N S(J) = -S(J)*S(J-1) 116 CONTINUE WTNM = WTS+WTF DO 117 I=ITSP,ITFM WTNM = WTNM+SINT(I) 117 CONTINUE YPS = CZR*WTNM*(S(2)-1.) C C BOUNDARY CONDITION AT R=RF C 118 GO TO (119,120,120,119,119,120),NBDCND 119 CR = CN(N) JRF = N GO TO 121 120 CR = CN(N+1) AN(N+1) = AN(N+1)+CN(N+1) JRF = N+1 121 WRS = AN(JRS+1)*R(JRS)**2/CN(JRS) WRF = CN(JRF-1)*R(JRF)**2/AN(JRF) WRZ = AN(JRS)/CZR JRSP = JRS+1 JRFM = JRF-1 MUNK = ITF-ITS+1 NUNK = JRF-JRS+1 DO 122 I=ITS,ITF BMH(I) = BM(I) 122 CONTINUE ISING = 0 GO TO (132,132,123,132,132,123),NBDCND 123 GO TO (132,132,124,132,132,124,132,124,124),MBDCND 124 IF (ELMBDA) 132,125,125 125 ISING = 1 SUM = WTS*WRS+WTS*WRF+WTF*WRS+WTF*WRF IF (ICTR) 126,127,126 126 SUM = SUM+WRZ 127 DO 129 J=JRSP,JRFM R2 = R(J)**2 DO 128 I=ITSP,ITFM SUM = SUM+R2*SINT(I) 128 CONTINUE 129 CONTINUE DO 130 J=JRSP,JRFM SUM = SUM+(WTS+WTF)*R(J)**2 130 CONTINUE DO 131 I=ITSP,ITFM SUM = SUM+(WRS+WRF)*SINT(I) 131 CONTINUE HNE = SUM 132 GO TO (133,133,133,133,134,134,133,133,134),MBDCND 133 BM(ITS) = BMH(ITS)+ELMBDA/SINT(ITS)**2 134 GO TO (135,135,135,135,135,135,136,136,136),MBDCND 135 BM(ITF) = BMH(ITF)+ELMBDA/SINT(ITF)**2 136 DO 137 I=ITSP,ITFM BM(I) = BMH(I)+ELMBDA/SINT(I)**2 137 CONTINUE GO TO (138,138,140,140,142,142,138,140,142),MBDCND 138 DO 139 J=JRS,JRF F(2,J) = F(2,J)-AT*F(1,J)/R(J)**2 139 CONTINUE GO TO 142 140 DO 141 J=JRS,JRF F(1,J) = F(1,J)+TDT*BDTS(J)*AT/R(J)**2 141 CONTINUE 142 GO TO (143,145,145,143,143,145,147,147,147),MBDCND 143 DO 144 J=JRS,JRF F(M,J) = F(M,J)-CT*F(M+1,J)/R(J)**2 144 CONTINUE GO TO 147 145 DO 146 J=JRS,JRF F(M+1,J) = F(M+1,J)-TDT*BDTF(J)*CT/R(J)**2 146 CONTINUE 147 GO TO (151,151,153,153,148,148),NBDCND 148 IF (MBDCND-3) 155,149,155 149 YHLD = F(ITS,1)-CZR/TDT*(SIN(TF)*BDTF(2)-SIN(TS)*BDTS(2)) DO 150 I=1,MP1 F(I,1) = YHLD 150 CONTINUE GO TO 155 151 RS2 = (RS+DR)**2 DO 152 I=ITS,ITF F(I,2) = F(I,2)-AR*F(I,1)/RS2 152 CONTINUE GO TO 155 153 DO 154 I=ITS,ITF F(I,1) = F(I,1)+TDR*BDRS(I)*AR/RS**2 154 CONTINUE 155 GO TO (156,158,158,156,156,158),NBDCND 156 RF2 = (RF-DR)**2 DO 157 I=ITS,ITF F(I,N) = F(I,N)-CR*F(I,N+1)/RF2 157 CONTINUE GO TO 160 158 DO 159 I=ITS,ITF F(I,N+1) = F(I,N+1)-TDR*BDRF(I)*CR/RF**2 159 CONTINUE 160 CONTINUE PERTRB = 0. IF (ISING) 161,170,161 161 SUM = WTS*WRS*F(ITS,JRS)+WTS*WRF*F(ITS,JRF)+WTF*WRS*F(ITF,JRS)+ 1 WTF*WRF*F(ITF,JRF) IF (ICTR) 162,163,162 162 SUM = SUM+WRZ*F(ITS,1) 163 DO 165 J=JRSP,JRFM R2 = R(J)**2 DO 164 I=ITSP,ITFM SUM = SUM+R2*SINT(I)*F(I,J) 164 CONTINUE 165 CONTINUE DO 166 J=JRSP,JRFM SUM = SUM+R(J)**2*(WTS*F(ITS,J)+WTF*F(ITF,J)) 166 CONTINUE DO 167 I=ITSP,ITFM SUM = SUM+SINT(I)*(WRS*F(I,JRS)+WRF*F(I,JRF)) 167 CONTINUE PERTRB = SUM/HNE DO 169 J=1,NP1 DO 168 I=1,MP1 F(I,J) = F(I,J)-PERTRB 168 CONTINUE 169 CONTINUE 170 DO 172 J=JRS,JRF RSQ = R(J)**2 DO 171 I=ITS,ITF F(I,J) = RSQ*F(I,J) 171 CONTINUE 172 CONTINUE IFLG = INTL 173 CALL BLKTRI (IFLG,NP,NUNK,AN(JRS),BN(JRS),CN(JRS),MP,MUNK, 1 AM(ITS),BM(ITS),CM(ITS),IDIMF,F(ITS,JRS),IERROR,W) IFLG = IFLG+1 IF (IFLG-1) 174,173,174 174 IF (NBDCND) 177,175,177 175 DO 176 I=1,MP1 F(I,JRF+1) = F(I,JRS) 176 CONTINUE 177 IF (MBDCND) 180,178,180 178 DO 179 J=1,NP1 F(ITF+1,J) = F(ITS,J) 179 CONTINUE 180 XP = 0. IF (ICTR) 181,188,181 181 IF (ISING) 186,182,186 182 SUM = WTS*F(ITS,2)+WTF*F(ITF,2) DO 183 I=ITSP,ITFM SUM = SUM+SINT(I)*F(I,2) 183 CONTINUE YPH = CZR*SUM XP = (F(ITS,1)-YPH)/YPS DO 185 J=JRS,JRF XPS = XP*S(J) DO 184 I=ITS,ITF F(I,J) = F(I,J)+XPS 184 CONTINUE 185 CONTINUE 186 DO 187 I=1,MP1 F(I,1) = XP 187 CONTINUE 188 RETURN END SUBROUTINE POIS (IFLG,NPEROD,N,MPEROD,M,A,B,C,IDIMY,Y,IERROR,W) POI02647 C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE POIS SOLVES THE LINEAR SYSTEM OF EQUATIONS C C A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J) C C + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J) C C FOR I = 1,2,...,M AND J = 1,2,...,N. C C THE INDICES I+1 AND I-1 ARE EVALUATED MODULO M, I.E., C X(0,J) = X(M,J) AND X(M+1,J) = X(1,J), AND X(I,0) MAY BE EQUAL TO C 0, X(I,2), OR X(I,N) AND X(I,N+1) MAY BE EQUAL TO 0, X(I,N-1), OR C X(I,1) DEPENDING ON AN INPUT PARAMETER. C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C IFLG C = 0 ON INITIAL ENTRY TO POIS OR IF N AND NPEROD ARE CHANGED C FROM PREVIOUS CALL. C = 1 IF N AND NPEROD ARE UNCHANGED FROM PREVIOUS CALL TO POIS. C C NOTE: A CALL WITH IFLG = 1 IS ABOUT 1 PERCENT FASTER THAN A C CALL WITH IFLG = 0 . C C NPEROD C INDICATES THE VALUES WHICH X(I,0) AND X(I,N+1) ARE ASSUMED TO C HAVE. C C = 0 IF X(I,0) = X(I,N) AND X(I,N+1) = X(I,1). C = 1 IF X(I,0) = X(I,N+1) = 0 . C = 2 IF X(I,0) = 0 AND X(I,N+1) = X(I,N-1). C = 3 IF X(I,0) = X(I,2) AND X(I,N+1) = X(I,N-1). C C N C THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. N IS DEPENDENT ON C NPEROD AND MUST HAVE THE FOLLOWING FORM: C C NPEROD C = 0 OR 2, THEN N = (2**P)(3**Q)(5**R) C = 1, THEN N = (2**P)(3**Q)(5**R)-1 C = 3, THEN N = (2**P)(3**Q)(5**R)+1 C C WHERE P, Q, AND R MAY BE ANY NON-NEGATIVE INTEGERS. N MUST C BE GREATER THAN 2 . C C MPEROD C = 0 IF A(1) AND C(M) ARE NOT ZERO C = 1 IF A(1) = C(M) = 0 C C M C THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. M MAY BE ANY INTEGER C GREATER THAN 1 . C C A,B,C C ONE-DIMENSIONAL ARRAYS OF LENGTH M WHICH SPECIFY THE C COEFFICIENTS IN THE LINEAR EQUATIONS GIVEN ABOVE. C C IDIMY C THE ROW (OR FIRST) DIMENSION OF THE TWO-DIMENSIONAL ARRAY Y AS C IT APPEARS IN THE PROGRAM CALLING POIS. THIS PARAMETER IS USED C TO SPECIFY THE VARIABLE DIMENSION OF Y. IDIMY MUST BE AT LEAST C M. C C Y C A TWO-DIMENSIONAL ARRAY WHICH SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE LINEAR SYSTEM OF EQUATIONS GIVEN ABOVE. Y MUST BE C DIMENSIONED AT LEAST M*N. C C W C A ONE-DIMENSIONAL ARRAY-WHICH MUST BE PROVIDED BY THE USER FOR C WORK SPACE. THE LENGTH OF W MUST BE AT LEAST 6N+5M. C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C Y C CONTAINS THE SOLUTION X. C C IERROR C AN ERROR FLAG WHICH INDICATES INVALID INPUT PARAMETERS EXCEPT C FOR NUMBER ZERO, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 M .LE. 2 . C = 2 N .LE. 2 OR NOT OF THE RIGHT FORM. C = 3 IDIMY .LT. M. C = 4 NPEROD .LT. 0 OR NPEROD .GT. 3 . C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C POIS WILL BE CALLED AGAIN WITH IFLAG = 1 . C C C EXTERNAL TRID ,TRIDP DIMENSION Y(IDIMY,1) DIMENSION W(1) ,B(1) ,A(1) ,C(1) IERROR = 0 IF (M .LE. 2) IERROR = 1 I = NCHECK(NPEROD-2*((2*NPEROD)/3)+N) IF (N.LE.2 .OR. I.GT.1) IERROR = 2 IF (IDIMY .LT. M) IERROR = 3 IF (NPEROD.LT.0 .OR. NPEROD.GT.3) IERROR = 4 IF (IERROR .NE. 0) GO TO 105 IWDIM1 = 6*N+1 IWDIM2 = IWDIM1+M IWDIM3 = IWDIM2+M IWDIM4 = IWDIM3+M IWDIM5 = IWDIM4+M DO 101 I=1,M A(I) = -A(I) C(I) = -C(I) K = IWDIM5+I-1 W(K) = -B(I)+2. 101 CONTINUE IF (MPEROD .EQ. 0) GO TO 102 CALL POISGN (NPEROD,N,IFLG,M,A,W(IWDIM5),C,IDIMY,Y,W(1), 1 W(IWDIM1),W(IWDIM2),W(IWDIM3),W(IWDIM4),TRID) GO TO 103 102 CALL POISGN (NPEROD,N,IFLG,M,A,W(IWDIM5),C,IDIMY,Y,W(1), 1 W(IWDIM1),W(IWDIM2),W(IWDIM3),W(IWDIM4),TRIDP) 103 DO 104 I=1,M A(I) = -A(I) C(I) = -C(I) 104 CONTINUE 105 RETURN END SUBROUTINE POISGN (NPEROD,N,IFLG,M,BA,BB,BC,IDIMQ,Q,TWOCOS,B,T,D, POI02801 1 W,TRI) DIMENSION Q(IDIMQ,1) DIMENSION BA(1) ,BB(1) ,BC(1) ,TWOCOS(1) , 1 B(1) ,T(1) ,D(1) ,W(1) IF (IFLG .NE. 0) GO TO 101 C C DO INITIALIZATION IF FIRST TIME THROUGH. C CALL POINIT (NPEROD,N,IEX2,IEX3,IEX5,N2PW,N2P3P,N2P3P5,KS3,KS5, 1 NPSTOP,TWOCOS) 101 MM1 = M-1 DO 104 J=1,N DO 102 I=1,M B(I) = -Q(I,J) 102 CONTINUE CALL TRI (1,1,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 103 I=1,M Q(I,J) = B(I) 103 CONTINUE 104 CONTINUE NP = 1 C C START REDUCTION FOR POWERS OF TWO C IF (IEX2 .EQ. 0) GO TO 112 L = IEX2 IF (IEX3+IEX5 .NE. 0) GO TO 105 IF (NPEROD .EQ. 1) L = L-1 IF (L .EQ. 0) GO TO 138 105 DO 111 KOUNT=1,L K = NP NP = 2*NP K2 = NP-1 K3 = NP K4 = 2*NP-1 JSTART = NP IF (NPEROD .EQ. 3) JSTART = 1 DO 110 J=JSTART,N,NP JM1 = J-K JP1 = J+K IF (J .EQ. 1) JM1 = JP1 IF (J .NE. N) GO TO 106 JP1 = JM1 IF (NPEROD .EQ. 0) JP1 = K 106 DO 107 I=1,M B(I) = Q(I,JM1)+Q(I,JP1) 107 CONTINUE CALL TRI (K,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 108 I=1,M T(I) = Q(I,J)+B(I) B(I) = T(I) 108 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 109 I=1,M Q(I,J) = T(I)+2.*B(I) 109 CONTINUE 110 CONTINUE 111 CONTINUE C C START REDUCTION FOR POWERS OF THREE C 112 IF (IEX3 .EQ. 0) GO TO 124 L = IEX3 IF (IEX5 .NE. 0) GO TO 113 IF (NPEROD .EQ. 1) L = L-1 IF (L .EQ. 0) GO TO 138 113 K2 = NP-1 DO 123 KOUNT=1,L K = NP NP = 3*NP K1 = K2+1 K2 = K2+K K3 = K2+1 K4 = K2+NP JSTART = NP IF (NPEROD .EQ. 3) JSTART = 1 DO 122 J=JSTART,N,NP IF (J .NE. 1) GO TO 114 JM1 = J+K JM2 = JM1+K GO TO 116 114 JM1 = J-K JM2 = JM1-K IF (J .NE. N) GO TO 116 IF (NPEROD .EQ. 0) GO TO 115 JP1 = JM1 JP2 = JM2 GO TO 117 115 JP1 = K JP2 = JP1+K GO TO 117 116 JP1 = J+K JP2 = JP1+K 117 DO 118 I=1,M B(I) = 2.*Q(I,J)+Q(I,JM2)+Q(I,JP2) 118 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 119 I=1,M T(I) = B(I)+Q(I,JM1)+Q(I,JP1) B(I) = T(I) 119 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 120 I=1,M Q(I,J) = Q(I,J)+B(I) B(I) = T(I) 120 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 121 I=1,M Q(I,J) = Q(I,J)+3.*B(I) 121 CONTINUE 122 CONTINUE 123 CONTINUE C C START REDUCTION FOR POWERS OF FIVE C 124 L = IEX5 IF (NPEROD .EQ. 1) L = L-1 IF (L .LE. 0) GO TO 138 K2 = (N2PW+N2P3P)/2-1 DO 137 KOUNT=1,L K = NP NP = 5*NP K1 = K2+1 K2 = K2+K K3 = K2+1 K4 = K2+NP JSTART = NP IF (NPEROD .EQ. 3) JSTART = 1 DO 136 J=JSTART,N,NP IF (J .NE. 1) GO TO 125 JM1 = J+K JM2 = JM1+K JM3 = JM2+K JM4 = JM3+K GO TO 127 125 JM1 = J-K JM2 = JM1-K JM3 = JM2-K JM4 = JM3-K IF (J .NE. N) GO TO 127 IF (NPEROD .EQ. 0) GO TO 126 JP1 = JM1 JP2 = JM2 JP3 = JM3 JP4 = JM4 GO TO 128 126 JP1 = K JP2 = JP1+K JP3 = JP2+K JP4 = JP3+K GO TO 128 127 JP1 = J+K JP2 = JP1+K JP3 = JP2+K JP4 = JP3+K 128 DO 129 I=1,M B(I) = 6.*Q(I,J)+4.*(Q(I,JM2)+Q(I,JP2))+Q(I,JM4)+Q(I,JP4) 129 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 130 I=1,M B(I) = B(I)+3.*(Q(I,JM1)+Q(I,JP1))+Q(I,JM3)+Q(I,JP3) 130 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 131 I=1,M T(I) = B(I) B(I) = 2.*Q(I,J)+Q(I,JM2)+Q(I,JP2)+B(I) 131 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 132 I=1,M B(I) = B(I)+Q(I,JM1)+Q(I,JP1) 132 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 133 I=1,M TEMP = B(I)+Q(I,J) B(I) = 4.*Q(I,J)+3.*(Q(I,JM2)+Q(I,JP2))+Q(I,JM4)+ 1 Q(I,JP4)-T(I) Q(I,J) = TEMP 133 CONTINUE CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 134 I=1,M B(I) = B(I)+2.*(Q(I,JM1)+Q(I,JP1))+Q(I,JM3)+Q(I,JP3) 134 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 135 I=1,M Q(I,J) = Q(I,J)+5.*B(I) 135 CONTINUE 136 CONTINUE 137 CONTINUE C C INITIAL PHASE OF BACK SUBSTITUTION C 138 IF (NPEROD.EQ.1 .OR. NPEROD.EQ.2) GO TO 147 IF (NPEROD .EQ. 0) GO TO 141 DO 139 I=1,M B(I) = 2.*Q(I,1) 139 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 140 I=1,M Q(I,N) = Q(I,N)+B(I) B(I) = 4.*Q(I,N) 140 CONTINUE CALL TRI (K4+1,K4+2*NP-1,2,M,MM1,BA,BB,BC,B,TWOCOS,D,W) GO TO 143 141 DO 142 I=1,M B(I) = 2.*Q(I,N) 142 CONTINUE CALL TRI (K4+1,K4+NP-1,2,M,MM1,BA,BB,BC,B,TWOCOS,D,W) 143 DO 144 I=1,M Q(I,N) = Q(I,N)+B(I) 144 CONTINUE IF (NPEROD .EQ. 0) GO TO 147 DO 145 I=1,M B(I) = 2.*Q(I,N) 145 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 146 I=1,M Q(I,1) = Q(I,1)+B(I) 146 CONTINUE C C REGULAR BACK SUBSTITUTION FOR POWERS OF FIVE C 147 CONTINUE IF (IEX5 .EQ. 0) GO TO 171 NP = N2P3P5 K8 = KS5 IF (NPEROD .EQ. 1) K3 = K3+NP/5 DO 170 KOUNT=1,IEX5 K = NP NP = NP/5 K4 = K3-1 K3 = K3-NP K1 = K8+1 K2 = K8+2*NP K5 = K2+1 K6 = K2+4*NP K7 = K6+1 K8 = K6+2*NP JSTART = NP IF (NPEROD .EQ. 3) JSTART = 1+NP DO 169 J=JSTART,N,K JM1 = J-NP JP1 = J+NP JP2 = JP1+NP JP3 = JP2+NP JP4 = JP3+NP IF (JM1 .NE. 0) GO TO 151 IF (NPEROD .EQ. 0) GO TO 149 DO 148 I=1,M B(I) = Q(I,JP1) 148 CONTINUE GO TO 153 149 DO 150 I=1,M B(I) = Q(I,JP1)+Q(I,N) 150 CONTINUE GO TO 153 151 DO 152 I=1,M B(I) = Q(I,JP1)+Q(I,JM1) 152 CONTINUE 153 CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 154 I=1,M Q(I,J) = Q(I,J)+B(I) B(I) = Q(I,JP1)+Q(I,JP3) 154 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 155 I=1,M Q(I,JP2) = Q(I,JP2)+B(I) 155 CONTINUE IF (JP4 .GT. N) GO TO 157 DO 156 I=1,M B(I) = 2.*Q(I,JP2)+Q(I,J)+Q(I,JP4) 156 CONTINUE GO TO 159 157 DO 158 I=1,M B(I) = 2.*Q(I,JP2)+Q(I,J) 158 CONTINUE 159 CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 160 I=1,M Q(I,JP2) = Q(I,JP2)+B(I) B(I) = Q(I,JP2)+Q(I,J) 160 CONTINUE CALL TRI (K5,K6,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 161 I=1,M Q(I,JP2) = Q(I,JP2)+B(I) B(I) = Q(I,J)+Q(I,JP2) 161 CONTINUE CALL TRI (K7,K8,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 162 I=1,M Q(I,J) = Q(I,J)+B(I) B(I) = Q(I,J)+Q(I,JP2) 162 CONTINUE CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 163 I=1,M Q(I,JP1) = Q(I,JP1)+B(I) 163 CONTINUE IF (JP4 .GT. N) GO TO 165 DO 164 I=1,M B(I) = Q(I,JP2)+Q(I,JP4) 164 CONTINUE GO TO 167 165 DO 166 I=1,M B(I) = Q(I,JP2) 166 CONTINUE 167 CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 168 I=1,M Q(I,JP3) = Q(I,JP3)+B(I) 168 CONTINUE 169 CONTINUE 170 CONTINUE C C REGULAR BACK SUBSTITUTION FOR POWERS OF THREE C 171 IF (IEX3 .EQ. 0) GO TO 191 NP = N2P3P K2 = KS3 IF (NPEROD.EQ.1 .AND. IEX5.EQ.0) K3 = K3+NP/3 DO 190 KOUNT=1,IEX3 K = NP NP = NP/3 K4 = K3-1 K3 = K3-NP K1 = K2+1 K2 = K2+2*NP JSTART = NP IF (NPEROD .EQ. 3) JSTART = NP+1 DO 189 J=JSTART,N,K JM1 = J-NP JP1 = J+NP JP2 = JP1+NP IF (JM1 .EQ. 0) GO TO 173 DO 172 I=1,M B(I) = Q(I,JP1)+Q(I,JM1) 172 CONTINUE GO TO 177 173 IF (NPEROD .EQ. 0) GO TO 175 DO 174 I=1,M B(I) = Q(I,JP1) 174 CONTINUE GO TO 177 175 DO 176 I=1,M B(I) = Q(I,JP1)+Q(I,N) 176 CONTINUE 177 CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 178 I=1,M Q(I,J) = Q(I,J)+B(I) 178 CONTINUE IF (JP2 .GT. N) GO TO 180 DO 179 I=1,M B(I) = Q(I,J)+Q(I,JP2) 179 CONTINUE GO TO 182 180 DO 181 I=1,M B(I) = Q(I,J) 181 CONTINUE 182 CALL TRI (K1,K2,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 183 I=1,M Q(I,J) = Q(I,J)+B(I) 183 CONTINUE IF (JP2 .GT. N) GO TO 185 DO 184 I=1,M B(I) = Q(I,J)+Q(I,JP2) 184 CONTINUE GO TO 187 185 DO 186 I=1,M B(I) = Q(I,J) 186 CONTINUE 187 CALL TRI (K3,K4,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 188 I=1,M Q(I,JP1) = Q(I,JP1)+B(I) 188 CONTINUE 189 CONTINUE 190 CONTINUE C C REGULAR BACK SUBSTITUTION FOR POWERS OF TWO C 191 IF (IEX2 .EQ. 0) GO TO 202 NP = N2PW DO 201 KOUNT=1,IEX2 K = NP NP = NP/2 K3 = K-1 JSTART = NP IF (NPEROD .EQ. 3) JSTART = 1+NP DO 200 J=JSTART,N,K JM1 = J-NP JP1 = J+NP IF (JM1 .NE. 0) GO TO 194 IF (JP1 .GT. N) GO TO 200 IF (NPEROD.EQ.1 .OR. NPEROD.EQ.2) GO TO 192 JM1 = N GO TO 196 192 DO 193 I=1,M B(I) = Q(I,JP1) 193 CONTINUE GO TO 198 194 IF (JP1 .LE. N) GO TO 196 DO 195 I=1,M B(I) = Q(I,JM1) 195 CONTINUE GO TO 198 196 DO 197 I=1,M B(I) = Q(I,JM1)+Q(I,JP1) 197 CONTINUE 198 CALL TRI (NP,K3,1,M,MM1,BA,BB,BC,B,TWOCOS,D,W) DO 199 I=1,M Q(I,J) = Q(I,J)+B(I) 199 CONTINUE 200 CONTINUE 201 CONTINUE 202 RETURN END SUBROUTINE POINIT (NPEROD,N,IEX2,IEX3,IEX5,N2PW,N2P3P,N2P3P5,KS3, POI03213 1 KS5,NPSTOP,TWOCOS) DIMENSION TWOCOS(1) C C PARAMETER NPSTOP IS NOT USED IN THIS SUBROUTINE. C C C MACHINE DEPENDENT CONSTANT C C PI=3.1415926535897932384626433832795028841971693993751058209749446 C PI = 3.14159265358979 C C COMPUTE EXPONENTS OF 2,3, AND 5 IN N. C NP = N IF (NPEROD .EQ. 1) NP = NP+1 IF (NPEROD .EQ. 3) NP = NP-1 IEX2 = 0 101 K = NP/2 IF (2*K .NE. NP) GO TO 102 IEX2 = IEX2+1 NP = K GO TO 101 102 IEX3 = 0 103 K = NP/3 IF (3*K .NE. NP) GO TO 104 IEX3 = IEX3+1 NP = K GO TO 103 104 IEX5 = 0 105 K = NP/5 IF (5*K .NE. NP) GO TO 106 IEX5 = IEX5+1 NP = K GO TO 105 106 CONTINUE N2PW = 2**IEX2 N2P3P = N2PW*(3**IEX3) N2P3P5 = N2P3P*(5**IEX5) C C COMPUTE NECESSARY VALUES OF 2*COS(X) C NP = 1 TWOCOS(1) = 0. K = 1 IF (IEX2 .EQ. 0) GO TO 110 L = IEX2 IF (IEX3+IEX5 .NE. 0) GO TO 107 IF (NPEROD .EQ. 1) L = L-1 IF (L .EQ. 0) GO TO 129 107 DO 109 KOUNT=1,L NP = 2*NP DO 108 I=1,NP J = K+I TWOCOS(J) = 2.*COS((FLOAT(I)-.5)*PI/FLOAT(NP)) 108 CONTINUE K = K+NP 109 CONTINUE 110 IF (IEX3 .EQ. 0) GO TO 114 L = IEX3 IF (IEX5 .NE. 0) GO TO 111 IF (NPEROD .EQ. 1) L = L-1 IF (L .EQ. 0) GO TO 117 111 DO 113 KOUNT=1,L NP = 3*NP DO 112 I=1,NP J = K+I TWOCOS(J) = 2.*COS((FLOAT(I)-.5)*PI/FLOAT(NP)) 112 CONTINUE K = K+NP 113 CONTINUE 114 L = IEX5 IF (NPEROD .EQ. 1) L = L-1 IF (L .LE. 0) GO TO 117 DO 116 KOUNT=1,L NP = 5*NP DO 115 I=1,NP J = K+I TWOCOS(J) = 2.*COS((FLOAT(I)-.5)*PI/FLOAT(NP)) 115 CONTINUE K = K+NP 116 CONTINUE 117 IF (NPEROD.EQ.1 .OR. NPEROD.EQ.2) GO TO 121 IF (NPEROD .EQ. 0) GO TO 119 NPT2 = 2*NP DO 118 I=1,NPT2 J = K+I TWOCOS(J) = 2.*COS(FLOAT(I)*PI/FLOAT(NP)) 118 CONTINUE K = K+NPT2 GO TO 121 119 DO 120 I=1,NP J = K+I TWOCOS(J) = 2.*COS(2.*FLOAT(I)*PI/FLOAT(NP)) 120 CONTINUE K = K+NP 121 NP = N2P3P5 IF (IEX5 .EQ. 0) GO TO 126 KS5 = K DO 125 KOUNT=1,IEX5 NP = NP/5 NPT2 = 2*NP DO 122 I=1,NPT2 J = K+I TWOCOS(J) = 2.*COS((FLOAT(I)-.5)*PI/FLOAT(NPT2)) 122 CONTINUE K = K+NPT2 DO 123 I=1,NP J = K+4*I TWOCOS(J-3) = 2.*COS((FLOAT(I)-.8)*PI/FLOAT(NP)) TWOCOS(J-2) = 2.*COS((FLOAT(I)-.6)*PI/FLOAT(NP)) TWOCOS(J-1) = 2.*COS((FLOAT(I)-.4)*PI/FLOAT(NP)) TWOCOS(J) = 2.*COS((FLOAT(I)-.2)*PI/FLOAT(NP)) 123 CONTINUE K = K+4*NP DO 124 I=1,NP J = K+2*I TWOCOS(J-1) = 2.*COS(FLOAT(3*I-2)*PI/FLOAT(3*NP)) TWOCOS(J) = 2.*COS(FLOAT(3*I-1)*PI/FLOAT(3*NP)) 124 CONTINUE K = K+2*NP 125 CONTINUE 126 IF (IEX3 .EQ. 0) GO TO 129 KS3 = K DO 128 KOUNT=1,IEX3 NP = NP/3 DO 127 I=1,NP J = K+2*I TWOCOS(J-1) = 2.*COS(FLOAT(3*I-2)*PI/FLOAT(3*NP)) TWOCOS(J) = 2.*COS(FLOAT(3*I-1)*PI/FLOAT(3*NP)) 127 CONTINUE K = K+2*NP 128 CONTINUE 129 RETURN END SUBROUTINE TRID (KSTART,KSTOP,ISING,M,MM1,A,B,C,Y,TWOCOS,D,W) POI03350 DIMENSION A(1) ,B(1) ,C(1) ,Y(1) , 1 TWOCOS(1) ,D(1) ,W(1) C C SUBROUTINE TO SOLVE TRIDIAGONAL SYSTEMS C C C PARAMETER W NOT USED IN THIS SUBROUTINE. C DO 103 K=KSTART,KSTOP X = -TWOCOS(K) D(1) = C(1)/(B(1)+X) Y(1) = Y(1)/(B(1)+X) DO 101 I=2,M Z = B(I)+X-A(I)*D(I-1) D(I) = C(I)/Z Y(I) = (Y(I)-A(I)*Y(I-1))/Z 101 CONTINUE DO 102 IP=1,MM1 I = M-IP Y(I) = Y(I)-D(I)*Y(I+1) 102 CONTINUE 103 CONTINUE IF (ISING .EQ. 1) RETURN D(1) = C(1)/(B(1)-2.) Y(1) = Y(1)/(B(1)-2.) DO 104 I=2,MM1 Z = B(I)-2.-A(I)*D(I-1) D(I) = C(I)/Z Y(I) = (Y(I)-A(I)*Y(I-1))/Z 104 CONTINUE Z = B(M)-2.-A(M)*D(M-1) X = Y(M)-A(M)*Y(M-1) IF (Z .NE. 0.) GO TO 105 Y(M) = 0. GO TO 106 105 Y(M) = X/Z 106 DO 107 IP=1,MM1 I = M-IP Y(I) = Y(I)-D(I)*Y(I+1) 107 CONTINUE RETURN END SUBROUTINE TRIDP (KSTART,KSTOP,ISING,M,MM1,A,B,C,Y,TWOCOS,D,W) POI03394 DIMENSION A(1) ,B(1) ,C(1) ,Y(1) , 1 TWOCOS(1) ,D(1) ,W(1) C C SUBROUTINE TO SOLVE PERIODIC TRIDIAGONAL SYSTEM C DO 103 K=KSTART,KSTOP X = -TWOCOS(K) D(1) = C(1)/(B(1)+X) W(1) = A(1)/(B(1)+X) Y(1) = Y(1)/(B(1)+X) BM = B(M) Z = C(M) DO 101 I=2,MM1 DEN = B(I)+X-A(I)*D(I-1) D(I) = C(I)/DEN W(I) = -A(I)*W(I-1)/DEN Y(I) = (Y(I)-A(I)*Y(I-1))/DEN Y(M) = Y(M)-Z*Y(I-1) BM = BM-Z*W(I-1) Z = -Z*D(I-1) 101 CONTINUE D(MM1) = D(MM1)+W(MM1) Z = A(M)+Z DEN = BM+X-Z*D(MM1) Y(M) = Y(M)-Z*Y(M-1) Y(M) = Y(M)/DEN Y(MM1) = Y(MM1)-D(MM1)*Y(M) DO 102 IP=2,MM1 I = M-IP Y(I) = Y(I)-D(I)*Y(I+1)-W(I)*Y(M) 102 CONTINUE 103 CONTINUE IF (ISING .EQ. 1) RETURN D(1) = C(1)/(B(1)-2.) W(1) = A(1)/(B(1)-2.) Y(1) = Y(1)/(B(1)-2.) BM = B(M) Z = C(M) DO 104 I=2,MM1 DEN = B(I)-2.-A(I)*D(I-1) D(I) = C(I)/DEN W(I) = -A(I)*W(I-1)/DEN Y(I) = (Y(I)-A(I)*Y(I-1))/DEN Y(M) = Y(M)-Z*Y(I-1) BM = BM-Z*W(I-1) Z = -Z*D(I-1) 104 CONTINUE D(MM1) = D(MM1)+W(MM1) Z = A(M)+Z DEN = BM-2.-Z*D(MM1) Y(M) = Y(M)-Z*Y(M-1) IF (DEN .NE. 0.) GO TO 105 Y(M) = 0. GO TO 106 105 Y(M) = Y(M)/DEN 106 Y(MM1) = Y(MM1)-D(MM1)*Y(M) DO 107 IP=2,MM1 I = M-IP Y(I) = Y(I)-D(I)*Y(I+1)-W(I)*Y(M) 107 CONTINUE RETURN END SUBROUTINE BLKTRI (IFLG,NP,N,AN,BN,CN,MP,M,AM,BM,CM,IDIMY,Y, POI03458 1 IERROR,W) C C C*********************************************************************** C C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 C C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN C C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS C C BY C C PAUL SWARZTRAUBER AND ROLAND SWEET C C TECHNICAL NOTE TN/IA-109 JULY 1975 C C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307 C C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION C C*********************************************************************** C C C C SUBROUTINE BLKTRI SOLVES A SYSTEM OF LINEAR EQUATIONS OF THE FORM C C AN(J)*X(I,J-1) + AM(I)*X(I-1,J) + (BN(J)+BM(I))*X(I,J) C C + CN(J)*X(I,J+1) + CM(I)*X(I+1,J) = Y(I,J) C C FOR I = 1,2,...,M AND J = 1,2,...,N. C C I+1 AND I-1 ARE EVALUATED MODULO M AND J+1 AND J-1 MODULO N, I.E., C C X(I,0) = X(I,N), X(I,N+1) = X(I,1), C X(0,J) = X(M,J), X(M+1,J) = X(1,J). C C THESE EQUATIONS USUALLY RESULT FROM THE DISCRETIZATION OF C SEPARABLE ELLIPTIC EQUATIONS. BOUNDARY CONDITIONS MAY BE C DIRICHLET, NEUMANN, OR PERIODIC. C C C * * * * * * * * * * ON INPUT * * * * * * * * * * C C IFLG C = 0 INITIALIZATION ONLY. CERTAIN QUANTITIES THAT DEPEND ON NP, C N, AN, BN, AND CN ARE COMPUTED AND C STORED IN THE WORK ARRAY W. C = 1 THE QUANTITIES THAT WERE COMPUTED IN THE INITIALIZATION ARE C USED TO OBTAIN THE SOLUTION X(I,J). C C NOTE A CALL WITH IFLG=0 TAKES APPROXIMATELY ONE HALF THE TIME C TIME AS A CALL WITH IFLG = 1 . HOWEVER, THE C INITIALIZATION DOES NOT HAVE TO BE REPEATED UNLESS NP, N, C AN, BN, OR CN CHANGE. C C NP C = 0 IF AN(1) AND CN(N) ARE NOT ZERO, WHICH CORRESPONDS TO C PERIODIC BOUNARY CONDITIONS. C = 1 IF AN(1) AND CN(N) ARE ZERO. C C N C THE NUMBER OF UNKNOWNS IN THE J-DIRECTION. IF NP = 1, N MUST BE C OF THE FORM 2**K-1 WHERE K IS AN INTEGER .GT. 1 . IF NP = 0, N C MUST BE OF THE FORM 2**K. (THE OPERATION COUNT OF THE ALGORITHM C IS PROPORTIONAL TO MN LOG2N AND, HENCE, N SHOULD BE SELECTED C LESS THAN OR EQUAL TO M.) C C AN,BN,CN C ONE-DIMENSIONAL ARRAYS OF LENGTH N THAT SPECIFY THE COEFFICIENTS C IN THE LINEAR EQUATIONS GIVEN ABOVE. C C MP C = 0 IF AM(1) AND CM(M) ARE NOT ZERO, WHICH CORRESPONDS TO C PERIODIC BOUNDARY CONDITIONS. C = 1 IF AM(1) = CM(M) = 0 . C C M C THE NUMBER OF UNKNOWNS IN THE I-DIRECTION. M MAY BE ANY INTEGER C GREATER THAN 4. C C AM,BM,CM C ONE-DIMENSIONAL ARRAYS OF LENGTH M THAT SPECIFY THE COEFFICIENTS C IN THE LINEAR EQUATIONS GIVEN ABOVE. C C IDIMY C THE ROW (OR FIRST) DIMENSION OF THE TWO-DIMENSIONAL ARRAY Y AS C IT APPEARS IN THE PROGRAM CALLING BLKTRI. THIS PARAMETER IS C USED TO SPECIFY THE VARIABLE DIMENSION OF Y. IDIMY MUST BE AT C LEAST M. C C Y C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE VALUES OF THE RIGHT C SIDE OF THE LINEAR SYSTEM OF EQUATIONS GIVEN ABOVE. Y MUST BE C DIMENSIONED AT LEAST M*N. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE PROVIDED BY THE USER FOR C WORK SPACE. IF NP = 0, THE LENGTH OF W MUST BE AT LEAST C (2NLOG2(N)+N+2+MAX(4N,6M)). C IF NP = 1, (2(N+1)(LOG2(N+1)-1)+2+MAX(2N,6M)). C C C * * * * * * * * * * ON OUTPUT * * * * * * * * * * C C Y C CONTAINS THE SOLUTION X. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT PARAMETERS. EXCEPT C FOR NUMBER ZERO, A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR. C = 1 M .LT. 5 . C = 2 N IS NOT THE FORM 2**K-1 WHEN NP = 1 . C = 3 N IS NOT OF THE FORM 2**K WHEN NP = 0 . C = 4 BLKTRI FAILED WHILE COMPUTING RESULTS THAT DEPEND ON THE C COEFFICIENT ARRAYS AN, BN, CN. CHECK THESE ARRAYS. C = 5 IDIMY.LT. M. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT BE DESTROYED IF C BLKTRI WILL BE CALLED AGAIN WITH IFLG = 1 . C C C DIMENSION AN(1) ,BN(1) ,CN(1) ,AM(1) , 1 BM(1) ,CM(1) ,Y(IDIMY,1) ,W(1) EXTERNAL PROD ,PRODP ,CPROD ,CPRODP COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NCMPLX ,IK ,IZ C C TEST M AND N FOR THE PROPER FORM C IERROR = 1 IF (M-5) 117,101,101 101 NH = N NPP = NP IF (NPP) 102,103,102 102 NH = NH+1 103 IK = 2 K = 0 104 IK = IK+IK K = K+1 IF (NH/IK) 105,105,104 105 IF (NPP) 106,108,106 106 IERROR = 2 IWCN = 2*(N+1)*(K-1)-N+3 IW1 = IWCN+N IWAH = IW1 IWBH = IWAH+N IF (K-2) 117,107,107 107 NCK = 2**K-1 IF (N-NCK) 117,110,117 108 IERROR = 3 IWCN = 2*K*N+3 IW1 = IWCN+N IWAH = IW1 IWBH = IWAH+N+N IF (K-2) 117,109,109 109 NCK = 2**K IF (N-NCK) 117,110,117 C C DIVIDE W INTO SEVERAL SUB WORKING ARRAYS C 110 IERROR = 5 IF (IDIMY-M) 117,111,111 111 IERROR = 0 IW2 = IW1+M IW3 = IW2+M IWD = IW3+M IWW = IWD+M IWU = IWW+M IF (IFLG) 114,112,114 112 W(IWCN) = CN(N) I = IWCN DO 113 J=2,N I = I+1 W(I) = CN(J-1) 113 CONTINUE C C SUBROUTINE COMP B COMPUTES THE ROOTS OF THE B POLYNOMIALS C CALL COMPB (N,IERROR,AN,BN,W(IWCN),W,W(IWAH),W(IWBH)) GO TO 117 114 IF (MP) 115,116,115 C C SUBROUTINE BLKTR1 SOLVES THE LINEAR SYSTEM C 115 CALL BLKTR1 (N,AN,BN,W(IWCN),M,AM,BM,CM,IDIMY,Y,W,W(IW1),W(IW2), 1 W(IW3),W(IWD),W(IWW),W(IWU),PROD,CPROD) GO TO 117 116 CALL BLKTR1 (N,AN,BN,W(IWCN),M,AM,BM,CM,IDIMY,Y,W,W(IW1),W(IW2), 1 W(IW3),W(IWD),W(IWW),W(IWU),PRODP,CPRODP) 117 CONTINUE RETURN END SUBROUTINE BLKTR1 (N,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,B,W1,W2,W3,WD, POI03659 1 WW,WU,PRDCT,CPRDCT) C C BLKTR1 SOLVES THE LINEAR SYSTEM C C B CONTAINS THE ROOTS OF ALL THE B POLYNOMIALS C W1,W2,W3,WD,WW,WU ARE ALL WORKING ARRAYS C PRDCT IS EITHER PRODP OR PROD DEPENDING ON WHETHER THE BOUNDARY C CONDITIONS IN THE M DIRECTION ARE PERIODIC OR NOT C CPRDCT IS EITHER CPRODP OR CPROD WHICH ARE THE COMPLEX VERSIONS C OF PRODP AND PROD. THESE ARE CALLED IN THE EVENT THAT SOME C OF THE ROOTS OF THE B SUB P POLYNOMIAL ARE COMPLEX C C DIMENSION AN(1) ,BN(1) ,CN(1) ,AM(1) , 1 BM(1) ,CM(1) ,B(1) ,W1(1) , 2 W2(1) ,W3(1) ,WD(1) ,WW(1) , 3 WU(1) ,Y(IDIMY,1) COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NCMPLX ,IK ,IZ NH = 2**K IH1 = NH+NH C C BEGIN REDUCTION PHASE C KDO = K IF (NPP) 101,102,101 101 KDO = K-1 102 DO 108 L=1,KDO IR = L-1 IZ = 2**IR ISGN = (-1)**IR MSGN = -ISGN IH2 = 2**(K-IR+1) LM = (IR-2)*IH1+IH2+IH2+1 LZ = (IR-1)*IH1+IH2+1 IIM = IZ-1 IIZ = IIM+IIM+1 JM1 = IIM+IIM+LM I1 = IZ+IZ CALL PRDCT (IIZ,B(LZ),IIM,B(LM),IIM,B(JM1),0,DUM,Y(1,IZ),W3,M, 1 AM,BM,CM,WD,WW,WU,ISGN) IF = 2**(K-IR-1) DO 107 JJ=1,IF I = JJ*I1 I6 = I I7 = I-IZ I9 = I+IZ J2 = JJ+JJ J4 = J2+J2 JM1 = (J4-2)*IIM+LM JP1 = J4*IIM+LM JP2 = J2*IIZ+LZ JP3 = (J4+2)*IIM+LM IF (JJ-IF) 105,103,103 103 IF (NPP) 107,104,107 104 JP1 = LM JP2 = LZ JP3 = IIM+IIM+LM I6 = 0 I9 = IZ 105 CALL PRDCT (IIM,B(JM1),0,DUM,0,DUM,IZ,AN(I7+1),W3,W1,M,AM, 1 BM,CM,WD,WW,WU,MSGN) CALL PRDCT (IIZ,B(JP2),IIM,B(JP1),IIM,B(JP3),0,DUM,Y(1,I9), 1 W3,M,AM,BM,CM,WD,WW,WU,ISGN) CALL PRDCT (IIM,B(JP1),0,DUM,0,DUM,IZ,CN(I6+1),W3,W2,M,AM, 1 BM,CM,WD,WW,WU,MSGN) DO 106 J=1,M Y(J,I) = W1(J)+W2(J)-Y(J,I) 106 CONTINUE 107 CONTINUE 108 CONTINUE IF (NPP) 112,109,112 109 J2 = 2*N*(K-1)+3 J1 = 2*N*(K-2)+5 IF (NCMPLX) 110,111,110 110 CALL CPRDCT (N,B(J2),N-1,B(J1),0,DUM,0,DUM,Y(1,N),Y(1,N),M,AM,BM, 1 CM,W1,W3,WW,MSGN) GO TO 112 111 CALL PRDCT (N,B(J2),N-1,B(J1),0,DUM,0,DUM,Y(1,N),Y(1,N),M,AM,BM, 1 CM,WD,WW,WU,MSGN) C C BEGIN BACK SUBSTITUTION PHASE C 112 DO 126 LL=1,K L = K-LL+1 IR = L-1 ISGN = (-1)**IR MSGN = -ISGN IZ = 2**IR IIM = IZ-1 IIZ = IIM+IIM+1 IH2 = 2**(K-IR+1) LM = (IR-2)*IH1+IH2+IH2+1 LZ = (IR-1)*IH1+IH2+1 IF = 2**(K-IR)-1 DO 125 JJ=1,IF,2 I = JJ*IZ I5 = I-IZ I6 = I+IZ I7 = I5 J2 = JJ+JJ JM1 = (J2-2)*IIM+LM JZ = (JJ-1)*IIZ+LZ JP1 = J2*IIM+LM IF (JJ-1) 113,113,117 113 IF (NPP) 115,114,115 114 I7 = N GO TO 117 115 DO 116 J=1,M W1(J) = 0. 116 CONTINUE GO TO 118 117 CALL PRDCT (IIM,B(JM1),0,DUM,0,DUM,IZ,AN(I5+1),Y(1,I7),W1, 1 M,AM,BM,CM,WD,WW,WU,MSGN) 118 IF (JJ-IF) 122,119,119 119 IF (NPP) 120,122,120 120 DO 121 J=1,M W2(J) = 0. 121 CONTINUE GO TO 123 122 CALL PRDCT (IIM,B(JP1),0,DUM,0,DUM,IZ,CN(I+1),Y(1,I6),W2,M, 1 AM,BM,CM,WD,WW,WU,MSGN) 123 DO 124 J=1,M W1(J) = Y(J,I)-W1(J)-W2(J) 124 CONTINUE CALL PRDCT (IIZ,B(JZ),IIM,B(JM1),IIM,B(JP1),0,DUM,W1, 1 Y(1,I),M,AM,BM,CM,WD,WW,WU,ISGN) 125 CONTINUE 126 CONTINUE RETURN END SUBROUTINE PROD (ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,W,U,IS)POI03792 C C PROD APPLIES A SEQUENCE OF MATRIX OPERATIONS TO THE VECTOR X AND C STORES THE RESULT IN Y C BD,BM1,BM2 ARE ARRAYS CONTAINING ROOTS OF CERTIAN B POLYNOMIALS C ND,NM1,NM2 ARE THE LENGTHS OF THE ARRAYS BD,BM1,BM2 RESPECTIVELY C AA ARRAY CONTAINING SCALAR MULTIPLIERS OF THE VECTOR X C NA IS THE LENGTH OF THE ARRAY AA C X,Y THE MATRIX OPERATIONS ARE APPLIED TO X AND THE RESULT IS Y C A,B,C ARE ARRAYS WHICH CONTAIN THE TRIDIAGONAL MATRIX C M IS THE ORDER OF THE MATRIX C D,W,U ARE WORKING ARRAYS C IS DETERMINES WHETHER OR NOT A CHANGE IN SIGN IS MADE C DIMENSION A(1) ,B(1) ,C(1) ,X(1) , 1 Y(1) ,D(1) ,W(1) ,BD(1) , 2 BM1(1) ,BM2(1) ,AA(1) ,U(1) IF (ND) 102,102,101 101 IF (IS) 104,104,102 102 DO 103 J=1,M W(J) = X(J) Y(J) = X(J) 103 CONTINUE GO TO 106 104 DO 105 J=1,M W(J) = -X(J) Y(J) = W(J) 105 CONTINUE 106 MM = M-1 ID = ND IBR = 0 M1 = NM1 M2 = NM2 IA = NA 107 IF (IA) 110,110,108 108 RT = AA(IA) IA = IA-1 C C SCALAR MULTIPLICATION C DO 109 J=1,M Y(J) = RT*W(J) 109 CONTINUE 110 IF (ID) 130,130,111 111 RT = BD(ID) ID = ID-1 IF (ID .EQ. 0) IBR = 1 C C BEGIN SOLUTION TO SYSTEM C D(M) = A(M)/(B(M)-RT) W(M) = Y(M)/(B(M)-RT) DO 112 J=2,MM K = M-J DEN = B(K+1)-RT-C(K+1)*D(K+2) D(K+1) = A(K+1)/DEN W(K+1) = (Y(K+1)-C(K+1)*W(K+2))/DEN 112 CONTINUE DEN = B(1)-RT-C(1)*D(2) W(1) = 1. IF (DEN) 113,114,113 113 W(1) = (Y(1)-C(1)*W(2))/DEN 114 DO 115 J=2,M W(J) = W(J)-D(J)*W(J-1) 115 CONTINUE IF (NA) 118,118,107 116 DO 117 J=1,M Y(J) = W(J) 117 CONTINUE IBR = 1 GO TO 107 118 IF (M1) 119,119,120 119 IF (M2) 116,116,125 120 IF (M2) 122,122,121 121 IF (ABS(BM1(M1))-ABS(BM2(M2))) 125,125,122 122 IF (IBR) 123,123,124 123 IF (ABS(BM1(M1)-BD(ID))-ABS(BM1(M1)-RT)) 116,124,124 124 RT = RT-BM1(M1) M1 = M1-1 GO TO 128 125 IF (IBR) 126,126,127 126 IF (ABS(BM2(M2)-BD(ID))-ABS(BM2(M2)-RT)) 116,127,127 127 RT = RT-BM2(M2) M2 = M2-1 128 DO 129 J=1,M Y(J) = Y(J)+RT*W(J) 129 CONTINUE GO TO 107 130 RETURN END SUBROUTINE PRODP (ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,U,W, POI03883 1 IS) C C PRODP APPLIES A SEQUENCE OF MATRIX OPERATIONS TO THE VECTOR X AND C STORES THE RESULT IN Y PERIODIC BOUNDARY CONDITIONS C C BD,BM1,BM2 ARE ARRAYS CONTAINING ROOTS OF CERTIAN B POLYNOMIALS C ND,NM1,NM2 ARE THE LENGTHS OF THE ARRAYS BD,BM1,BM2 RESPECTIVELY C AA ARRAY CONTAINING SCALAR MULTIPLIERS OF THE VECTOR X C NA IS THE LENGTH OF THE ARRAY AA C X,Y THE MATRIX OPERATIONS ARE APPLIED TO X AND THE RESULT IS Y C A,B,C ARE ARRAYS WHICH CONTAIN THE TRIDIAGONAL MATRIX C M IS THE ORDER OF THE MATRIX C D,U,W ARE WORKING ARRAYS C IS DETERMINES WHETHER OR NOT A CHANGE IN SIGN IS MADE C DIMENSION A(1) ,B(1) ,C(1) ,X(1) , 1 Y(1) ,D(1) ,U(1) ,BD(1) , 2 BM1(1) ,BM2(1) ,AA(1) ,W(1) IF (ND) 102,102,101 101 IF (IS) 104,104,102 102 DO 103 J=1,M Y(J) = X(J) W(J) = X(J) 103 CONTINUE GO TO 106 104 DO 105 J=1,M Y(J) = -X(J) W(J) = Y(J) 105 CONTINUE 106 MM = M-1 MM2 = M-2 ID = ND IBR = 0 M1 = NM1 M2 = NM2 IA = NA 107 IF (IA) 110,110,108 108 RT = AA(IA) IA = IA-1 DO 109 J=1,M Y(J) = RT*W(J) 109 CONTINUE 110 IF (ID) 131,131,111 111 RT = BD(ID) ID = ID-1 IF (ID .EQ. 0) IBR = 1 C C BEGIN SOLUTION TO SYSTEM C BH = B(M)-RT YM = Y(M) DEN = B(1)-RT D(1) = C(1)/DEN U(1) = A(1)/DEN W(1) = Y(1)/DEN V = C(M) DO 112 J=2,MM2 DEN = B(J)-RT-A(J)*D(J-1) D(J) = C(J)/DEN U(J) = -A(J)*U(J-1)/DEN W(J) = (Y(J)-A(J)*W(J-1))/DEN BH = BH-V*U(J-1) YM = YM-V*W(J-1) V = -V*D(J-1) 112 CONTINUE DEN = B(M-1)-RT-A(M-1)*D(M-2) D(M-1) = (C(M-1)-A(M-1)*U(M-2))/DEN W(M-1) = (Y(M-1)-A(M-1)*W(M-2))/DEN AM = A(M)-V*D(M-2) BH = BH-V*U(M-2) YM = YM-V*W(M-2) DEN = BH-AM*D(M-1) IF (DEN) 113,114,113 113 W(M) = (YM-AM*W(M-1))/DEN GO TO 115 114 W(M) = 1. 115 W(M-1) = W(M-1)-D(M-1)*W(M) DO 116 J=2,MM K = M-J W(K) = W(K)-D(K)*W(K+1)-U(K)*W(M) 116 CONTINUE IF (NA) 119,119,107 117 DO 118 J=1,M Y(J) = W(J) 118 CONTINUE IBR = 1 GO TO 107 119 IF (M1) 120,120,121 120 IF (M2) 117,117,126 121 IF (M2) 123,123,122 122 IF (ABS(BM1(M1))-ABS(BM2(M2))) 126,126,123 123 IF (IBR) 124,124,125 124 IF (ABS(BM1(M1)-BD(ID))-ABS(BM1(M1)-RT)) 117,125,125 125 RT = RT-BM1(M1) M1 = M1-1 GO TO 129 126 IF (IBR) 127,127,128 127 IF (ABS(BM2(M2)-BD(ID))-ABS(BM2(M2)-RT)) 117,128,128 128 RT = RT-BM2(M2) M2 = M2-1 129 DO 130 J=1,M Y(J) = Y(J)+RT*W(J) 130 CONTINUE GO TO 107 131 RETURN END SUBROUTINE CPROD (ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,YY,M,A,B,C,D,W,Y, POI03991 1 ISGN) C C PROD APPLIES A SEQUENCE OF MATRIX OPERATIONS TO THE VECTOR X AND C STORES THE RESULT IN YY (COMPLEX CASE) C AA ARRAY CONTAINING SCALAR MULTIPLIERS OF THE VECTOR X C ND,NM1,NM2 ARE THE LENGTHS OF THE ARRAYS BD,BM1,BM2 RESPECTIVELY C BD,BM1,BM2 ARE ARRAYS CONTAINING ROOTS OF CERTIAN B POLYNOMIALS C NA IS THE LENGTH OF THE ARRAY AA C X,YY THE MATRIX OPERATIONS ARE APPLIED TO X AND THE RESULT IS YY C A,B,C ARE ARRAYS WHICH CONTAIN THE TRIDIAGONAL MATRIX C M IS THE ORDER OF THE MATRIX C D,W,Y ARE WORKING ARRAYS C ISGN DETERMINES WHETHER OR NOT A CHANGE IN SIGN IS MADE C COMPLEX Y ,D ,W ,BD , 1 CRT ,DEN ,Y1 ,Y2 DIMENSION A(1) ,B(1) ,C(1) ,X(1) , 1 Y(1) ,D(1) ,W(1) ,BD(1) , 2 BM1(1) ,BM2(1) ,AA(1) ,YY(1) IF (ND) 102,102,101 101 IF (ISGN) 104,104,102 102 DO 103 J=1,M Y(J) = CMPLX(X(J),0.) 103 CONTINUE GO TO 106 104 DO 105 J=1,M Y(J) = CMPLX(-X(J),0.) 105 CONTINUE 106 MM = M-1 ID = ND M1 = NM1 M2 = NM2 IA = NA 107 IFLG = 0 IF (ID) 114,114,108 108 CRT = BD(ID) ID = ID-1 IFLG = 1 C C BEGIN SOLUTION TO SYSTEM C D(M) = A(M)/(B(M)-CRT) W(M) = Y(M)/(B(M)-CRT) DO 109 J=2,MM K = M-J DEN = B(K+1)-CRT-C(K+1)*D(K+2) D(K+1) = A(K+1)/DEN W(K+1) = (Y(K+1)-C(K+1)*W(K+2))/DEN 109 CONTINUE DEN = B(1)-CRT-C(1)*D(2) IF (CABS(DEN)) 110,111,110 110 Y(1) = (Y(1)-C(1)*W(2))/DEN GO TO 112 111 Y(1) = (1.,0.) 112 DO 113 J=2,M Y(J) = W(J)-D(J)*Y(J-1) 113 CONTINUE 114 IF (M1) 115,115,117 115 IF (M2) 126,126,116 116 RT = BM2(M2) M2 = M2-1 GO TO 122 117 IF (M2) 118,118,119 118 RT = BM1(M1) M1 = M1-1 GO TO 122 119 IF (ABS(BM1(M1))-ABS(BM2(M2))) 121,121,120 120 RT = BM1(M1) M1 = M1-1 GO TO 122 121 RT = BM2(M2) M2 = M2-1 122 Y1 = (B(1)-RT)*Y(1)+C(1)*Y(2) IF (MM-2) 125,123,123 C C MATRIX MULTIPLICATION C 123 DO 124 J=2,MM Y2 = A(J)*Y(J-1)+(B(J)-RT)*Y(J)+C(J)*Y(J+1) Y(J-1) = Y1 Y1 = Y2 124 CONTINUE 125 Y(M) = A(M)*Y(M-1)+(B(M)-RT)*Y(M) Y(M-1) = Y1 IFLG = 1 GO TO 107 126 IF (IA) 129,129,127 127 RT = AA(IA) IA = IA-1 IFLG = 1 C C SCALAR MULTIPLICATION C DO 128 J=1,M Y(J) = RT*Y(J) 128 CONTINUE 129 IF (IFLG) 130,130,107 130 DO 131 J=1,M YY(J) = REAL(Y(J)) 131 CONTINUE RETURN END SUBROUTINE CPRODP (ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,YY,M,A,B,C,D,U, POI04095 1 Y,ISGN) C C PRODP APPLIES A SEQUENCE OF MATRIX OPERATIONS TO THE VECTOR X AND C STORES THE RESULT IN YY PERIODIC BOUNDARY CONDITIONS C AND COMPLEX CASE C C BD,BM1,BM2 ARE ARRAYS CONTAINING ROOTS OF CERTIAN B POLYNOMIALS C ND,NM1,NM2 ARE THE LENGTHS OF THE ARRAYS BD,BM1,BM2 RESPECTIVELY C AA ARRAY CONTAINING SCALAR MULTIPLIERS OF THE VECTOR X C NA IS THE LENGTH OF THE ARRAY AA C X,YY THE MATRIX OPERATIONS ARE APPLIED TO X AND THE RESULT IS YY C A,B,C ARE ARRAYS WHICH CONTAIN THE TRIDIAGONAL MATRIX C M IS THE ORDER OF THE MATRIX C D,U,Y ARE WORKING ARRAYS C ISGN DETERMINES WHETHER OR NOT A CHANGE IN SIGN IS MADE C COMPLEX Y ,D ,U ,V , 1 DEN ,BH ,YM ,AM , 2 Y1 ,Y2 ,YH ,BD , 3 CRT DIMENSION A(1) ,B(1) ,C(1) ,X(1) , 1 Y(1) ,D(1) ,U(1) ,BD(1) , 2 BM1(1) ,BM2(1) ,AA(1) ,YY(1) IF (ND) 102,102,101 101 IF (ISGN) 104,104,102 102 DO 103 J=1,M Y(J) = CMPLX(X(J),0.) 103 CONTINUE GO TO 106 104 DO 105 J=1,M Y(J) = CMPLX(-X(J),0.) 105 CONTINUE 106 MM = M-1 MM2 = M-2 ID = ND M1 = NM1 M2 = NM2 IA = NA 107 IFLG = 0 IF (ID) 114,114,108 108 CRT = BD(ID) ID = ID-1 IFLG = 1 C C BEGIN SOLUTION TO SYSTEM C BH = B(M)-CRT YM = Y(M) DEN = B(1)-CRT D(1) = C(1)/DEN U(1) = A(1)/DEN Y(1) = Y(1)/DEN V = CMPLX(C(M),0.) DO 109 J=2,MM2 DEN = B(J)-CRT-A(J)*D(J-1) D(J) = C(J)/DEN U(J) = -A(J)*U(J-1)/DEN Y(J) = (Y(J)-A(J)*Y(J-1))/DEN BH = BH-V*U(J-1) YM = YM-V*Y(J-1) V = -V*D(J-1) 109 CONTINUE DEN = B(M-1)-CRT-A(M-1)*D(M-2) D(M-1) = (C(M-1)-A(M-1)*U(M-2))/DEN Y(M-1) = (Y(M-1)-A(M-1)*Y(M-2))/DEN AM = A(M)-V*D(M-2) BH = BH-V*U(M-2) YM = YM-V*Y(M-2) DEN = BH-AM*D(M-1) IF (CABS(DEN)) 110,111,110 110 Y(M) = (YM-AM*Y(M-1))/DEN GO TO 112 111 Y(M) = (1.,0.) 112 Y(M-1) = Y(M-1)-D(M-1)*Y(M) DO 113 J=2,MM K = M-J Y(K) = Y(K)-D(K)*Y(K+1)-U(K)*Y(M) 113 CONTINUE 114 IF (M1) 115,115,117 115 IF (M2) 126,126,116 116 RT = BM2(M2) M2 = M2-1 GO TO 122 117 IF (M2) 118,118,119 118 RT = BM1(M1) M1 = M1-1 GO TO 122 119 IF (ABS(BM1(M1))-ABS(BM2(M2))) 121,121,120 120 RT = BM1(M1) M1 = M1-1 GO TO 122 121 RT = BM2(M2) M2 = M2-1 C C MATRIX MULTIPLICATION C 122 YH = Y(1) Y1 = (B(1)-RT)*Y(1)+C(1)*Y(2)+A(1)*Y(M) IF (MM-2) 125,123,123 123 DO 124 J=2,MM Y2 = A(J)*Y(J-1)+(B(J)-RT)*Y(J)+C(J)*Y(J+1) Y(J-1) = Y1 Y1 = Y2 124 CONTINUE 125 Y(M) = A(M)*Y(M-1)+(B(M)-RT)*Y(M)+C(M)*YH Y(M-1) = Y1 IFLG = 1 GO TO 107 126 IF (IA) 129,129,127 127 RT = AA(IA) IA = IA-1 IFLG = 1 C C SCALAR MULTIPLICATION C DO 128 J=1,M Y(J) = RT*Y(J) 128 CONTINUE 129 IF (IFLG) 130,130,107 130 DO 131 J=1,M YY(J) = REAL(Y(J)) 131 CONTINUE RETURN END SUBROUTINE PPADD (N,IERROR,A,C,CBP,BP,BH,BN) POI04221 C C PPADD COMPUTES THE ROOTS OF THE B SUB P POLYNOMIAL C THIS ROUTINE IS CALLED AT THE LAST STEP OF THE PREPROCESSING PHASE C IN THE CASE OF PERIODIC BOUNDARY CONDITIONS C C N IS THE ORDER OF THE BH AND BP POLYNOMIALS C BP IS WHERE THE ROOTS OF THE B SUB P POLYNOMIAL ARE STORED C CBP IS THE SAME AS BP EXCEPT TYPE COMPLEX C BH IS USED TO TEMPORARILY STORE THE ROOTS OF THE B HAT POLYNOMIAL C WHICH ENTERS THROUGH BP C BN IS TEMPORARY STORAGE USED TO INDICATE THE TYPE OF ROOT IN BP C WHETHER REAL, DOUBLE OR COMPLEX C COMPLEX CF ,CX ,FSG ,HSG , 1 DD ,F ,FP ,FPP , 2 CDIS ,R1 ,R2 ,R3 , 3 CBP DIMENSION A(1) ,C(1) ,BP(1) ,BH(1) , 1 BN(1) ,CBP(1) COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NCMPLX ,IK ,IZ EXTERNAL PSGF ,PPSPF ,PPSGF SCNV = SQRT(CNV) IZ = N IZM2 = IZ-2 DO 101 J=1,IZ BH(J) = BP(J) 101 CONTINUE XL = BH(1) DB = BH(3)-BH(1) 102 XL = XL-DB IF (PSGF(XL,IZ,C,A,BH)) 102,103,103 103 SGN = -1. CBP(1) = CMPLX(BSRH(XL,BH(1),IZ,C,A,BH,PSGF,SGN),0.) XR = BH(IZ) DB = BH(IZ)-BH(IZ-2) 104 XR = XR+DB IF (PSGF(XR,IZ,C,A,BH)) 104,105,105 105 SGN = 1. CBP(IZ) = CMPLX(BSRH(BH(IZ),XR,IZ,C,A,BH,PSGF,SGN),0.) DO 121 IG=2,IZM2,2 XL = BH(IG) XR = BH(IG+1) SGN = -1. XM = BSRH(XL,XR,IZ,C,A,BH,PPSPF,SGN) PSG = PSGF(XM,IZ,C,A,BH) IF (PSG-EPS) 107,106,106 C C CASE OF A REAL ZERO C 106 SGN = 1. CBP(IG) = CMPLX(BSRH(BH(IG),XM,IZ,C,A,BH,PSGF,SGN),0.) SGN = -1. CBP(IG+1) = CMPLX(BSRH(XM,BH(IG+1),IZ,C,A,BH,PSGF,SGN),0.) BN(IG) = 0. BN(IG+1) = 0. GO TO 121 C C CASE OF A COMPLEX ZERO C 107 IT = 0 ICV = 0 CX = CMPLX(XM,0.) 108 FSG = (1.,0.) HSG = (1.,0.) FP = (0.,0.) FPP = (0.,0.) DO 109 J=1,IZ DD = 1./(CX-BH(J)) FSG = FSG*A(J)*DD HSG = HSG*C(J)*DD FP = FP+DD FPP = FPP-DD*DD 109 CONTINUE F = (1.,0.)-FSG-HSG I3 = 0 IF (CABS(FP)) 111,111,110 110 I3 = 1 R3 = -F/FP 111 I2 = 0 IF (CABS(FPP)) 117,117,112 112 I2 = 1 CDIS = CSQRT(FP**2-2.*F*FPP) R1 = CDIS-FP R2 = -FP-CDIS IF (CABS(R1)-CABS(R2)) 114,114,113 113 R1 = R1/FPP GO TO 115 114 R1 = R2/FPP 115 R2 = 2.*F/FPP/R1 IF (CABS(R2) .LT. CABS(R1)) R1 = R2 IF (I3) 118,118,116 116 IF (CABS(R3) .LT. CABS(R1)) R1 = R3 GO TO 118 117 R1 = R3 118 CX = CX+R1 IT = IT+1 IF (IT .GT. 50) GO TO 124 IF (CABS(R1) .GT. SCNV) GO TO 108 IF (ICV) 119,119,120 119 ICV = 1 GO TO 108 120 CBP(IG) = CX CBP(IG+1) = CONJG(CX) 121 CONTINUE NCMPLX = 1 DO 122 J=2,IZ IF (AIMAG(CBP(J))) 125,122,125 122 CONTINUE NCMPLX = 0 DO 123 J=2,IZ BP(J) = REAL(CBP(J)) 123 CONTINUE GO TO 125 124 IERROR = 4 125 CONTINUE RETURN END FUNCTION PSGF (X,IZ,C,A,BH) POI04341 DIMENSION A(1) ,C(1) ,BH(1) FSG = 1. HSG = 1. DO 101 J=1,IZ DD = 1./(X-BH(J)) FSG = FSG*A(J)*DD HSG = HSG*C(J)*DD 101 CONTINUE PSGF = 1.-FSG-HSG RETURN END FUNCTION BSRH (XLL,XRR,IZ,C,A,BH,F,SGN) POI04354 DIMENSION A(1) ,C(1) ,BH(1) COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NCMPLX ,IK ,IZZ XL = XLL XR = XRR DX = .5*ABS(XR-XL) 101 X = .5*(XL+XR) IF (SGN*F(X,IZ,C,A,BH)) 103,105,102 102 XR = X GO TO 104 103 XL = X 104 DX = .5*DX IF (DX-CNV) 105,105,101 105 BSRH = .5*(XL+XR) RETURN END FUNCTION PPSGF (X,IZ,C,A,BH) POI04372 DIMENSION A(1) ,C(1) ,BH(1) SUM = 0. DO 101 J=1,IZ SUM = SUM-1./(X-BH(J))**2 101 CONTINUE PPSGF = SUM RETURN END FUNCTION PPSPF (X,IZ,C,A,BH) POI04382 DIMENSION A(1) ,C(1) ,BH(1) SUM = 0. DO 101 J=1,IZ SUM = SUM+1./(X-BH(J)) 101 CONTINUE PPSPF = SUM RETURN END SUBROUTINE COMPB (N,IERROR,AN,BN,CN,B,AH,BH) POI04392 C C COMPB COMPUTES THE ROOTS OF THE B POLYNOMIALS USING THE EISPACK C SUBROUTINE IMTQL1. IERROR IS SET TO 4 IF IF EITHER IMTQL1 FAILS C OR IF A(J+1)*C(J) IS LESS THAN ZERO FOR SOME J. C AH,BH ARE TEMPORARY WORK ARRAYS C DIMENSION AN(1) ,BN(1) ,CN(1) ,B(1) , 1 AH(1) ,BH(1) COMMON /CBLKT/ NPP ,K ,EPS ,CNV , 1 NCMPLX ,IK ,IZ EPS = APXEPS(DUM) BNORM = ABS(BN(1)) DO 101 J=2,N BNORM = AMAX1(BNORM,ABS(BN(J))) 101 CONTINUE CNV = EPS*BNORM DO 107 J=1,N IF (J-1) 103,102,103 102 IF (NPP) 107,103,107 103 ARG = AN(J)*CN(J) IF (ARG) 126,104,104 104 CHLD = SQRT(ARG) IF (CN(J)) 106,105,105 105 B(J) = CHLD GO TO 107 106 B(J) = -CHLD 107 CONTINUE IF = 2**K LH = IF LH1 = IF+1 I1 = 1 DO 119 L=1,K I1 = I1+I1 NN = I1+I1-1 IFL = IF-I1+1 DO 118 I=1,IFL,I1 IF (I-IFL) 113,108,108 108 IF (NPP) 109,110,109 109 LH = LH+NN GO TO 118 110 LS = 0 DO 111 J=I,IF LS = LS+1 BH(LS) = BN(J) AH(LS) = B(J) 111 CONTINUE I1M = I1-1 DO 112 J=1,I1M LS = LS+1 BH(LS) = BN(J) AH(LS) = B(J) 112 CONTINUE GO TO 115 113 LS = 0 JFL = I+NN-1 DO 114 J=I,JFL LS = LS+1 BH(LS) = BN(J) AH(LS) = B(J) 114 CONTINUE 115 CALL TQLRAT (NN,BH,AH,IERROR) IF (IERROR) 126,116,126 116 DO 117 J=1,NN LH = LH+1 B(LH) = -BH(J) 117 CONTINUE 118 CONTINUE LH1 = LH+1 119 CONTINUE DO 120 J=1,N B(J) = -BN(J) 120 CONTINUE IF (NPP) 125,121,125 121 J2 = 2*N*(K-1)+4 LH = J2 J1 = J2-N-N+1 N2M2 = J2+N+N-4 122 D1 = ABS(B(J1)-B(J2-1)) D2 = ABS(B(J1)-B(J2)) D3 = ABS(B(J1)-B(J2+1)) IF ((D2 .LT. D1) .AND. (D2 .LT. D3)) GO TO 123 B(LH) = B(J2) J2 = J2+1 LH = LH+1 IF (J2-N2M2) 122,122,124 123 J2 = J2+1 J1 = J1+1 IF (J2-N2M2) 122,122,124 124 B(LH) = B(N2M2+1) J1 = 2*N*(K-1)+3 J2 = J1+N+N+N J3 = J2+N CALL PPADD (N,IERROR,AN,CN,B(J1),B(J1),B(J2),B(J3)) 125 RETURN 126 IERROR = 4 RETURN END SUBROUTINE TQLRAT (N,D,E2,IERR) POI04491 C INTEGER I ,J ,L ,M , 1 N ,II ,L1 ,MML , 2 IERR REAL D(N) ,E2(N) REAL B ,C ,F ,G , 1 H ,P ,R ,S , 2 MACHEP C C REAL SQRT,ABS,SIGN C COMMON /CBLKT/ NPP ,K ,MACHEP ,CNV , 1 NCMPLX ,IK ,IZ C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT- C C N IS THE ORDER OF THE MATRIX, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E2 CONTAINS THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES, C C E2 HAS BEEN DESTROYED, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** C IERR = 0 IF (N .EQ. 1) GO TO 114 C DO 101 I=2,N E2(I-1) = E2(I)*E2(I) 101 CONTINUE C F = 0.0 B = 0.0 E2(N) = 0.0 C DO 112 L=1,N J = 0 H = MACHEP*(ABS(D(L))+SQRT(E2(L))) IF (B .GT. H) GO TO 102 B = H C = B*B C C ********** LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ********** C 102 DO 103 M=L,N IF (E2(M) .LE. C) GO TO 104 C C ********** E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP ********** C 103 CONTINUE C 104 IF (M .EQ. L) GO TO 108 105 IF (J .EQ. 30) GO TO 113 J = J+1 C C ********** FORM SHIFT ********** C L1 = L+1 S = SQRT(E2(L)) G = D(L) P = (D(L1)-G)/(2.0*S) R = SQRT(P*P+1.0) D(L) = S/(P+SIGN(R,P)) H = G-D(L) C DO 106 I=L1,N D(I) = D(I)-H 106 CONTINUE C F = F+H C C ********** RATIONAL QL TRANSFORMATION ********** C G = D(M) IF (G .EQ. 0.0) G = B H = G S = 0.0 MML = M-L C C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** C DO 107 II=1,MML I = M-II P = G*H R = P+E2(I) E2(I+1) = S*R S = E2(I)/R D(I+1) = H+S*(H+D(I)) G = D(I)-E2(I)/G IF (G .EQ. 0.0) G = B H = G*P/R 107 CONTINUE C E2(L) = S*G D(L) = H C C ********** GUARD AGAINST UNDERFLOWED H ********** C IF (H .EQ. 0.0) GO TO 108 IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 108 E2(L) = H*E2(L) IF (E2(L) .NE. 0.0) GO TO 105 108 P = D(L)+F C C ********** ORDER EIGENVALUES ********** C IF (L .EQ. 1) GO TO 110 C C ********** FOR I=L STEP -1 UNTIL 2 DO -- ********** C ABSP = ABS(P) DO 109 II=2,L I = L+2-II IF (ABSP .GE. ABS(D(I-1))) GO TO 111 D(I) = D(I-1) 109 CONTINUE C 110 I = 1 111 D(I) = P 112 CONTINUE C GO TO 114 C C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** C 113 IERR = L 114 RETURN C C ********** LAST CARD OF TQLRAT ********** C END FUNCTION NCHECK (N) POI04654 C C THIS SUBPROGRAM CHECKS THE FORM OF N. C IF N = 2**P, THEN NCHECK = 0. C IF N = (2**P)(3**Q)(5**Q), THEN NCHECK = 1. C OTHERWISE, NCHECK = 2. C NP = N 101 K = NP/2 IF (2*K .NE. NP) GO TO 102 NP = K GO TO 101 102 IF (NP .NE. 1) GO TO 103 NCHECK = 0 RETURN 103 K = NP/3 IF (3*K .NE. NP) GO TO 104 NP = K GO TO 103 104 K = NP/5 IF (5*K .NE. NP) GO TO 105 NP = K GO TO 104 105 IF (NP .NE. 1) GO TO 106 NCHECK = 1 RETURN 106 NCHECK = 2 RETURN END FUNCTION APXEPS (DUM) POI04684 COMMON /VALUE/ V EPS = 1. 101 EPS = EPS/10. CALL STORE (EPS+1.) IF (V-1.) 102,102,101 102 APXEPS = 100.*EPS RETURN END SUBROUTINE STORE (X) POI04694 COMMON /VALUE/ V V = X RETURN END C PROGRAM XAMPLE POI04700 C POI04701 C***********************************************************************POI04702 C POI04703 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI04704 C POI04705 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI04706 C POI04707 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI04708 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI04709 C POI04710 C BY POI04711 C POI04712 C PAUL SWARZTRAUBER AND ROLAND SWEET POI04713 C POI04714 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI04715 C POI04716 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI04717 C POI04718 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI04719 C POI04720 C***********************************************************************POI04721 C POI04722 C POI04723 C PROGRAM TO ILLUSTRATE THE USE OF PWSCRT. POI04724 C POI04725 DIMENSION F(120,100) ,BDB(81) ,W(1294) ,X(101) , POI04726 1 Y(81) POI04727 C POI04728 C FROM DIMENSION STATEMENT WE GET VALUE OF IDIMF. ALSO NOTE THAT W POI04729 C IS DIMENSIONED 6*(N+1) + 8*(M+1). POI04730 C POI04731 IDIMF = 120 POI04732 INTL = 0 POI04733 A = 0. POI04734 B = 2. POI04735 M = 100 POI04736 MBDCND = 2 POI04737 C = -1. POI04738 D = 3. POI04739 N = 80 POI04740 NBDCND = 0 POI04741 ELMBDA = -4. POI04742 C POI04743 C AUXILIARY QUANTITIES. POI04744 C POI04745 PI = 3.14159265358979 POI04746 PIBY2 = PI/2. POI04747 PISQ = PI**2 POI04748 MP1 = M+1 POI04749 NP1 = N+1 POI04750 C POI04751 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING POI04752 C BOUNDARY DATA AND THE RIGHT SIDE OF THE HELMHOLTZ EQUATION. POI04753 C POI04754 DO 101 I=1,MP1 POI04755 X(I) = FLOAT(I-1)/50. POI04756 101 CONTINUE POI04757 DO 102 J=1,NP1 POI04758 Y(J) = -1.+FLOAT(J-1)/20. POI04759 102 CONTINUE POI04760 C POI04761 C GENERATE BOUNDARY DATA. POI04762 C POI04763 DO 103 J=1,NP1 POI04764 BDB(J) = 4.*COS((Y(J)+1.)*PIBY2) POI04765 103 CONTINUE POI04766 C POI04767 C BDA, BDC, AND BDD ARE DUMMY VARIABLES. POI04768 C POI04769 DO 104 J=1,NP1 POI04770 F(1,J) = 0. POI04771 104 CONTINUE POI04772 C POI04773 C GENERATE RIGHT SIDE OF EQUATION. POI04774 C POI04775 DO 106 I=2,MP1 POI04776 DO 105 J=1,NP1 POI04777 F(I,J) = (2.-(4.+PISQ/4.)*X(I)**2)*COS((Y(J)+1.)*PIBY2) POI04778 105 CONTINUE POI04779 106 CONTINUE POI04780 CALL PWSCRT (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, POI04781 1 ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI04782 C POI04783 C COMPUTE DISCRETIZATION ERROR. THE EXACT SOLUTION IS POI04784 C U(X,Y) = X**2*COS((Y+1)*PIBY2) POI04785 C POI04786 ERR = 0. POI04787 DO 108 I=1,MP1 POI04788 DO 107 J=1,NP1 POI04789 Z = ABS(F(I,J)-X(I)**2*COS((Y(J)+1.)*PIBY2)) POI04790 IF (Z .GT. ERR) ERR = Z POI04791 107 CONTINUE POI04792 108 CONTINUE POI04793 PRINT 1001 , IERROR,ERR POI04794 C POI04795 C POI04796 C POI04797 1001 FORMAT (///9H IERROR =,I2,10X,22HDISCRETIZATION ERROR =,E12.5) POI04798 C POI04799 END POI04800 C PROGRAM XAMPLE POI04802 C POI04803 C***********************************************************************POI04804 C POI04805 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI04806 C POI04807 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI04808 C POI04809 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI04810 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI04811 C POI04812 C BY POI04813 C POI04814 C PAUL SWARZTRAUBER AND ROLAND SWEET POI04815 C POI04816 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI04817 C POI04818 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI04819 C POI04820 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI04821 C POI04822 C***********************************************************************POI04823 C POI04824 C POI04825 C PROGRAM TO ILLUSTRATE THE USE OF PWSPLR. POI04826 C POI04827 DIMENSION F(100,50) ,BDC(51) ,BDD(51) ,W(702) , POI04828 1 R(51) ,THETA(49) POI04829 C POI04830 C FROM DIMENSION STATEMENT WE GET VALUE OF IDIMF. ALSO NOTE THAT W POI04831 C IS DIMENSIONED 6*(N+1) + 8*(M+1). POI04832 C POI04833 IDIMF = 100 POI04834 INTL = 0 POI04835 A = 0. POI04836 B = 1. POI04837 M = 50 POI04838 MBDCND = 5 POI04839 C = 0. POI04840 PI = 3.14159265358979 POI04841 D = PI/2. POI04842 N = 48 POI04843 NBDCND = 3 POI04844 ELMBDA = 0. POI04845 C POI04846 C AUXILIARY QUANTITIES. POI04847 C POI04848 MP1 = M+1 POI04849 NP1 = N+1 POI04850 C POI04851 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING POI04852 C BOUNDARY DATA AND THE RIGHT SIDE OF THE POISSON EQUATION. POI04853 C POI04854 DO 101 I=1,MP1 POI04855 R(I) = FLOAT(I-1)/50. POI04856 101 CONTINUE POI04857 DO 102 J=1,NP1 POI04858 THETA(J) = FLOAT(J-1)*PI/96. POI04859 102 CONTINUE POI04860 C POI04861 C GENERATE BOUNDARY DATA. POI04862 C POI04863 DO 103 I=1,MP1 POI04864 BDC(I) = 0. POI04865 BDD(I) = 0. POI04866 103 CONTINUE POI04867 C POI04868 C BDA AND BDB ARE DUMMY VARIABLES. POI04869 C POI04870 DO 104 J=1,NP1 POI04871 F(MP1,J) = 1.-COS(4.*THETA(J)) POI04872 104 CONTINUE POI04873 C POI04874 C GENERATE RIGHT SIDE OF EQUATION. POI04875 C POI04876 DO 106 I=1,M POI04877 DO 105 J=1,NP1 POI04878 F(I,J) = 16.*R(I)**2 POI04879 105 CONTINUE POI04880 106 CONTINUE POI04881 CALL PWSPLR (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, POI04882 1 ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI04883 C POI04884 C COMPUTE DISCRETIZATION ERROR. THE EXACT SOLUTION IS POI04885 C U(R,THETA) = R**4*(1 - COS(4*THETA)) POI04886 C POI04887 ERR = 0. POI04888 DO 108 I=1,MP1 POI04889 DO 107 J=1,NP1 POI04890 Z = ABS(F(I,J)-R(I)**4*(1.-COS(4.*THETA(J)))) POI04891 IF (Z .GT. ERR) ERR = Z POI04892 107 CONTINUE POI04893 108 CONTINUE POI04894 PRINT 1001 , IERROR,ERR POI04895 STOP POI04896 C POI04897 C POI04898 C POI04899 1001 FORMAT (///9H IERROR =,I2,10X,22HDISCRETIZATION ERROR =,E12.5) POI04900 C POI04901 END POI04902 C PROGRAM XAMPLE POI04904 C POI04905 C***********************************************************************POI04906 C POI04907 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI04908 C POI04909 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI04910 C POI04911 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI04912 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI04913 C POI04914 C BY POI04915 C POI04916 C PAUL SWARZTRAUBER AND ROLAND SWEET POI04917 C POI04918 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI04919 C POI04920 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI04921 C POI04922 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI04923 C POI04924 C***********************************************************************POI04925 C POI04926 C POI04927 C PROGRAM TO ILLUSTRATE THE USE OF PWSCYL. POI04928 C POI04929 DIMENSION F(75,105) ,BDA(101) ,BDB(101) ,BDC(51) , POI04930 1 BDD(51) ,W(1014) ,R(51) ,Z(101) POI04931 C POI04932 C FROM DIMENSION STATEMENT WE GET VALUE OF IDIMF. ALSO NOTE THAT W POI04933 C IS DIMENSIONED 6*(N+1) + 8*(M+1). POI04934 C POI04935 IDIMF = 75 POI04936 INTL = 0 POI04937 A = 0. POI04938 B = 1. POI04939 M = 50 POI04940 MBDCND = 6 POI04941 C = 0. POI04942 D = 1. POI04943 N = 100 POI04944 NBDCND = 3 POI04945 ELMBDA = 0. POI04946 C POI04947 C AUXILIARY QUANTITIES. POI04948 C POI04949 MP1 = M+1 POI04950 NP1 = N+1 POI04951 C POI04952 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING POI04953 C BOUNDARY DATA AND THE RIGHT SIDE OF THE POISSON EQUATION. POI04954 C POI04955 DO 101 I=1,MP1 POI04956 R(I) = FLOAT(I-1)/50. POI04957 101 CONTINUE POI04958 DO 102 J=1,NP1 POI04959 Z(J) = FLOAT(J-1)/100. POI04960 102 CONTINUE POI04961 C POI04962 C GENERATE BOUNDARY DATA. POI04963 C POI04964 DO 103 J=1,NP1 POI04965 BDB(J) = 4.*Z(J)**4 POI04966 103 CONTINUE POI04967 DO 104 I=1,MP1 POI04968 BDC(I) = 0. POI04969 BDD(I) = 4.*R(I)**4 POI04970 104 CONTINUE POI04971 C POI04972 C BDA IS A DUMMY VARIABLE. POI04973 C POI04974 C POI04975 C GENERATE RIGHT SIDE OF EQUATION. POI04976 C POI04977 DO 106 I=1,MP1 POI04978 DO 105 J=1,NP1 POI04979 F(I,J) = 4.*R(I)**2*Z(J)**2*(4.*Z(J)**2+3.*R(I)**2) POI04980 105 CONTINUE POI04981 106 CONTINUE POI04982 CALL PWSCYL (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC,BDD, POI04983 1 ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI04984 C POI04985 C COMPUTE DISCRETIZATION ERROR BY MINIMIZING OVER ALL A THE FUNCTIONPOI04986 C NORM(F(I,J) - A*1 - U(R(I),Z(J))). THE EXACT SOLUTION IS POI04987 C U(R,Z) = (R*Z)**4 + ARBITRARY CONSTANT. POI04988 C POI04989 X = 0. POI04990 DO 108 I=1,MP1 POI04991 DO 107 J=1,NP1 POI04992 X = X+F(I,J)-(R(I)*Z(J))**4 POI04993 107 CONTINUE POI04994 108 CONTINUE POI04995 X = X/FLOAT(NP1*MP1) POI04996 DO 110 I=1,MP1 POI04997 DO 109 J=1,NP1 POI04998 F(I,J) = F(I,J)-X POI04999 109 CONTINUE POI05000 110 CONTINUE POI05001 ERR = 0. POI05002 DO 112 I=1,MP1 POI05003 DO 111 J=1,NP1 POI05004 X = ABS(F(I,J)-(R(I)*Z(J))**4) POI05005 IF (X .GT. ERR) ERR = X POI05006 111 CONTINUE POI05007 112 CONTINUE POI05008 PRINT 1001 , PERTRB,IERROR,ERR POI05009 STOP POI05010 C POI05011 C POI05012 C POI05013 1001 FORMAT (///9H PERTRB =,E12.5,10X,8HIERROR =,I2,10X, POI05014 1 22HDISCRETIZATION ERROR =,E12.5) POI05015 C POI05016 END POI05017 C PROGRAM XAMPLE POI05019 C POI05020 C***********************************************************************POI05021 C POI05022 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI05023 C POI05024 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI05025 C POI05026 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI05027 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI05028 C POI05029 C BY POI05030 C POI05031 C PAUL SWARZTRAUBER AND ROLAND SWEET POI05032 C POI05033 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI05034 C POI05035 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI05036 C POI05037 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI05038 C POI05039 C***********************************************************************POI05040 C POI05041 C POI05042 C PROGRAM TO ILLUSTRATE THE USE OF PWSSSP POI05043 C POI05044 DIMENSION F(19,73) ,BDTF(73) ,SINT(19) ,SINP(73) , POI05045 1 W(647) POI05046 C POI05047 C THE VALUE OF IDIMF IS THE FIRST DIMENSION OF F. W IS DIMENSIONED POI05048 C 11*(M+1)+6X(N+1)=647 SINCE M=18 AND N=72 POI05049 C POI05050 PI = 3.141592653589793 POI05051 INTL = 0 POI05052 TS = 0 POI05053 TF = PI/2. POI05054 M = 18 POI05055 MBDCND = 6 POI05056 PS = 0 POI05057 PF = PI+PI POI05058 N = 72 POI05059 NBDCND = 0 POI05060 ELMBDA = 0. POI05061 IDIMF = 19 POI05062 C POI05063 C GENERATE SINES FOR USE IN SUBSEQUENT COMPUTATIONS POI05064 C POI05065 DTHETA = TF/FLOAT(M) POI05066 MP1 = M+1 POI05067 DO 101 I=1,MP1 POI05068 SINT(I) = SIN(FLOAT(I-1)*DTHETA) POI05069 101 CONTINUE POI05070 DPHI = (PI+PI)/FLOAT(N) POI05071 NP1 = N+1 POI05072 DO 102 J=1,NP1 POI05073 SINP(J) = SIN(FLOAT(J-1)*DPHI) POI05074 102 CONTINUE POI05075 C POI05076 C COMPUTE RIGHT SIDE OF EQUATION AND STORE IN F POI05077 C POI05078 DO 104 J=1,NP1 POI05079 DO 103 I=1,MP1 POI05080 F(I,J) = 2.-6.*(SINT(I)*SINP(J))**2 POI05081 103 CONTINUE POI05082 104 CONTINUE POI05083 C POI05084 C STORE DERIVATIVE DATA AT THE EQUATOR POI05085 C POI05086 DO 105 J=1,NP1 POI05087 BDTF(J) = 0. POI05088 105 CONTINUE POI05089 C POI05090 CALL PWSSSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS, POI05091 1 BDPF,ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI05092 C POI05093 C COMPUTE DISCRETIZATION ERROR. SINCE PROBLEM IS SINGULAR, THE POI05094 C SOLUTION MUST BE NORMALIZED. POI05095 C POI05096 ERR = 0 POI05097 DO 107 J=1,NP1 POI05098 DO 106 I=1,MP1 POI05099 Z = ABS(F(I,J)-(SINT(I)*SINP(J))**2-F(1,1)) POI05100 IF (Z .GT. ERR) ERR = Z POI05101 106 CONTINUE POI05102 107 CONTINUE POI05103 C POI05104 PRINT 1001 , IERROR,ERR,PERTRB POI05105 STOP POI05106 C POI05107 C POI05108 C POI05109 1001 FORMAT (///9H IERROR= ,I3,10X,22H DISCRETIZATION ERR = ,E12.4, POI05110 1 14H PERTURBATION=,E12.4) POI05111 C POI05112 END POI05113 C PROGRAM XAMPLE POI05115 C POI05116 C***********************************************************************POI05117 C POI05118 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI05119 C POI05120 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI05121 C POI05122 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI05123 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI05124 C POI05125 C BY POI05126 C POI05127 C PAUL SWARZTRAUBER AND ROLAND SWEET POI05128 C POI05129 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI05130 C POI05131 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI05132 C POI05133 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI05134 C POI05135 C***********************************************************************POI05136 C POI05137 C POI05138 C PROGRAM TO ILLUSTRATE THE USE OF PWSCSP POI05139 C POI05140 DIMENSION F(48,33) ,BDTF(33) ,W(902) ,R(33) , POI05141 1 THETA(48) POI05142 C POI05143 C THE VALUE OF IDIMF IS THE FIRST DIMENSION OF F. SINCE N=31,L=32 POI05144 C THEREFORE K=5 AND W IS DIMENSIONED 2(L+1)(K-1)+6(M+N)+MAX(4N,6M) POI05145 C +14 = 902 POI05146 C POI05147 PI = 3.141592653589793 POI05148 INTL = 0 POI05149 TS = 0. POI05150 TF = PI/2. POI05151 M = 36 POI05152 MBDCND = 6 POI05153 RS = 0. POI05154 RF = 1. POI05155 N = 32 POI05156 NBDCND = 5 POI05157 ELMBDA = 0. POI05158 IDIMF = 48 POI05159 C POI05160 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING THE POI05161 C BOUNDARY DATA AND THE RIGHT SIDE OF THE EQUATION. POI05162 C POI05163 MP1 = M+1 POI05164 DTHETA = TF/FLOAT(M) POI05165 DO 101 I=1,MP1 POI05166 THETA(I) = FLOAT(I-1)*DTHETA POI05167 101 CONTINUE POI05168 NP1 = N+1 POI05169 DR = 1./FLOAT(N) POI05170 DO 102 J=1,NP1 POI05171 R(J) = FLOAT(J-1)*DR POI05172 102 CONTINUE POI05173 C POI05174 C GENERATE NORMAL DERIVATIVE DATA AT EQUATOR POI05175 C POI05176 DO 103 J=1,NP1 POI05177 BDTF(J) = 0. POI05178 103 CONTINUE POI05179 C POI05180 C COMPUTE BOUNDARY DATA ON THE SURFACE OF THE SPHERE POI05181 C POI05182 DO 104 I=1,MP1 POI05183 F(I,N+1) = COS(THETA(I))**4 POI05184 104 CONTINUE POI05185 C POI05186 C COMPUTE RIGHT SIDE OF EQUATION POI05187 C POI05188 DO 106 I=1,MP1 POI05189 CI4 = 12.*COS(THETA(I))**2 POI05190 DO 105 J=1,N POI05191 F(I,J) = CI4*R(J)**2 POI05192 105 CONTINUE POI05193 106 CONTINUE POI05194 C POI05195 CALL PWSCSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND,BDRS, POI05196 1 BDRF,ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI05197 C POI05198 C COMPUTE DISCRETIZATION ERROR POI05199 C POI05200 ERR = 0. POI05201 DO 108 I=1,MP1 POI05202 CI4 = COS(THETA(I))**4 POI05203 DO 107 J=1,N POI05204 Z = ABS(F(I,J)-CI4*R(J)**4) POI05205 IF (Z .GT. ERR) ERR = Z POI05206 107 CONTINUE POI05207 108 CONTINUE POI05208 PRINT 1001 , IERROR,ERR POI05209 C POI05210 C THE FOLLOWING PROGRAM ILLUSTRATES THE USE OF PWSCSP TO SOLVE POI05211 C A THREE DIMENSIONAL PROBLEM WHICH HAS LONGITUDNAL DEPENDENCE POI05212 C POI05213 MBDCND = 2 POI05214 NBDCND = 1 POI05215 DPHI = PI/72. POI05216 ELMBDA = -2.*(1.-COS(DPHI))/DPHI**2 POI05217 C POI05218 C COMPUTE BOUNDARY DATA ON THE SURFACE OF THE SPHERE POI05219 C POI05220 DO 109 I=1,MP1 POI05221 F(I,N+1) = SIN(THETA(I)) POI05222 109 CONTINUE POI05223 C POI05224 C COMPUTE RIGHT SIDE OF THE EQUATION POI05225 C POI05226 DO 111 J=1,N POI05227 DO 110 I=1,MP1 POI05228 F(I,J) = 0. POI05229 110 CONTINUE POI05230 111 CONTINUE POI05231 C POI05232 CALL PWSCSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND,BDRS, POI05233 1 BDRF,ELMBDA,F,IDIMF,PERTRB,IERROR,W) POI05234 C POI05235 C COMPUTE DISCRETIZATION ERROR (FOURIER COEFFICIENTS) POI05236 C POI05237 ERR = 0 POI05238 DO 113 I=1,MP1 POI05239 SI = SIN(THETA(I)) POI05240 DO 112 J=1,NP1 POI05241 Z = ABS(F(I,J)-R(J)*SI) POI05242 IF (Z .GT. ERR) ERR = Z POI05243 112 CONTINUE POI05244 113 CONTINUE POI05245 C POI05246 PRINT 1001 , IERROR,ERR POI05247 STOP POI05248 C POI05249 C POI05250 C POI05251 1001 FORMAT (///9H IERROR =,I3,10X,22HDISCRETIZATION ERROR =,E12.4) POI05252 C POI05253 END POI05254 C PROGRAM XAMPLE POI05256 C POI05257 C***********************************************************************POI05258 C POI05259 C VERSION 3 OCTOBER 1978 POI05260 C POI05261 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI05262 C POI05263 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI05264 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI05265 C POI05266 C BY POI05267 C POI05268 C PAUL SWARZTRAUBER AND ROLAND SWEET POI05269 C POI05270 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI05271 C POI05272 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI05273 C POI05274 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI05275 C POI05276 C***********************************************************************POI05277 C POI05278 C POI05279 C PROGRAM TO ILLUSTRATE THE USE OF SUBROUTINE POIS. POI05280 C POI05281 DIMENSION F(25,130) ,A(20) ,B(20) ,C(20) , POI05282 1 W(820) ,X(20) ,Y(120) POI05283 C POI05284 C FROM DIMENSION STATEMENT WE GET VALUE OF IDIMY. ALSO NOTE THAT W POI05285 C IS DIMENSIONED 6*(N+1) + 5*(M+1). POI05286 C POI05287 IDIMY = 25 POI05288 IFLG = 0 POI05289 MPEROD = 1 POI05290 M = 20 POI05291 DELTAX = 1./FLOAT(M) POI05292 NPEROD = 0 POI05293 N = 120 POI05294 PI = 3.14159265358979 POI05295 DELTAY = 2.*PI/FLOAT(N) POI05296 C POI05297 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING POI05298 C COEFFICIENTS AND RIGHT SIDE OF EQUATION. POI05299 C POI05300 DO 101 I=1,M POI05301 X(I) = FLOAT(I-1)*DELTAX POI05302 101 CONTINUE POI05303 DO 102 J=1,120 POI05304 Y(J) = -PI+FLOAT(J-1)*DELTAY POI05305 102 CONTINUE POI05306 C POI05307 C GENERATE COEFFICIENTS. POI05308 C POI05309 S = (DELTAY/DELTAX)**2 POI05310 T = S*DELTAX POI05311 A(1) = 0. POI05312 B(1) = -2.*S POI05313 C(1) = 2.*S POI05314 DO 103 I=2,M POI05315 A(I) = (1.+X(I))**2*S+(1.+X(I))*T POI05316 C(I) = (1.+X(I))**2*S-(1.+X(I))*T POI05317 B(I) = -2.*(1.+X(I))**2*S POI05318 103 CONTINUE POI05319 C(M) = 0. POI05320 C POI05321 C GENERATE RIGHT SIDE OF EQUATION FOR I = 1 SHOWING INTRODUCTION OF POI05322 C BOUNDARY DATA. POI05323 C POI05324 DYSQ = DELTAY**2 POI05325 DO 104 J=1,N POI05326 F(1,J) = DYSQ*(11.+8./DELTAX)*SIN(Y(J)) POI05327 104 CONTINUE POI05328 C POI05329 C GENERATE RIGHT SIDE. POI05330 C POI05331 MM1 = M-1 POI05332 DO 106 I=2,MM1 POI05333 DO 105 J=1,N POI05334 F(I,J) = DYSQ*3.*(1.+X(I))**4*SIN(Y(J)) POI05335 105 CONTINUE POI05336 106 CONTINUE POI05337 C POI05338 C GENERATE RIGHT SIDE FOR I = M SHOWING INTRODUCTION OF POI05339 C BOUNDARY DATA. POI05340 C POI05341 DO 107 J=1,N POI05342 F(M,J) = DYSQ*(3.*(1.+X(M))**4-16.*((1.+X(M))/DELTAX)**2+ POI05343 1 16.*(1.+X(M))/DELTAX)*SIN(Y(J))POI05344 107 CONTINUE POI05345 CALL POIS (IFLG,NPEROD,N,MPEROD,M,A,B,C,IDIMY,F,IERROR,W) POI05346 C POI05347 C COMPUTE DISCRETIATION ERROR. THE EXACT SOLUTION IS POI05348 C U(X,Y) = (1+X)**4*SIN(Y) POI05349 C POI05350 ERR = 0. POI05351 DO 109 I=1,M POI05352 DO 108 J=1,N POI05353 Z = ABS(F(I,J)-(1.+X(I))**4*SIN(Y(J))) POI05354 IF (Z .GT. ERR) ERR = Z POI05355 108 CONTINUE POI05356 109 CONTINUE POI05357 PRINT 1001 , IERROR,ERR POI05358 STOP POI05359 C POI05360 C POI05361 C POI05362 1001 FORMAT (///9H IERROR =,I3,10X,22HDISCRETIZATION ERROR =,E12.4) POI05363 C POI05364 END POI05365 C PROGRAM XAMPLE POI05367 C POI05368 C***********************************************************************POI05369 C POI05370 C VERSION 2 OCTOBER 1976 INCLUDING ERRATA OCTOBER 1976 POI05371 C POI05372 C DOCUMENTATION FOR THIS PROGRAM IS GIVEN IN POI05373 C POI05374 C EFFICIENT FORTRAN SUBPROGRAMS FOR THE SOLUTION OF POI05375 C ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS POI05376 C POI05377 C BY POI05378 C POI05379 C PAUL SWARZTRAUBER AND ROLAND SWEET POI05380 C POI05381 C TECHNICAL NOTE TN/IA-109 JULY 1975 POI05382 C POI05383 C NATIONAL CENTER FOR ATMOSPHERIC RESEARCH BOULDER,COLORADO 80307POI05384 C POI05385 C WHICH IS SPONSORED BY THE NATIONAL SCIENCE FOUNDATION POI05386 C POI05387 C***********************************************************************POI05388 C POI05389 C POI05390 C PROGRAM TO ILLUSTRATE THE USE OF BLKTRI POI05391 C POI05392 DIMENSION Y(75,105) ,AM(75) ,BM(75) ,CM(75) , POI05393 1 AN(105) ,BN(105) ,CN(105) ,W(942) , POI05394 2 S(75) ,T(105) POI05395 C POI05396 C THE VALUE OF IDIMY IS THE FIRST DIMENSION OF Y. SINCE NP=1 W IS POI05397 C DIMENSIONED 2*(N+1)*(LOG2(N+1)-1)+2+MAX(2*N,6*M)=942 POI05398 C POI05399 IFLG = 0 POI05400 NP = 1 POI05401 N = 63 POI05402 MP = 1 POI05403 M = 50 POI05404 IDIMY = 75 POI05405 C POI05406 C GENERATE AND STORE GRID POINTS FOR THE PURPOSE OF COMPUTING THE POI05407 C COEFFICIENTS AND THE ARRAY Y. POI05408 C POI05409 DELTAS = 1./FLOAT(M+1) POI05410 DO 101 I=1,M POI05411 S(I) = FLOAT(I)*DELTAS POI05412 101 CONTINUE POI05413 DELTAT = 1./FLOAT(N+1) POI05414 DO 102 J=1,N POI05415 T(J) = FLOAT(J)*DELTAT POI05416 102 CONTINUE POI05417 C POI05418 C COMPUTE THE COEFFICIENTS AM,BM,CM CORRESPONDING TO THE S DIRECTIONPOI05419 C POI05420 HDS = DELTAS/2. POI05421 TDS = DELTAS+DELTAS POI05422 DO 103 I=1,M POI05423 TEMP1 = 1./(S(I)*TDS) POI05424 TEMP2 = 1./((S(I)-HDS)*TDS) POI05425 TEMP3 = 1./((S(I)+HDS)*TDS) POI05426 AM(I) = TEMP1*TEMP2 POI05427 CM(I) = TEMP1*TEMP3 POI05428 BM(I) = -(AM(I)+CM(I)) POI05429 103 CONTINUE POI05430 C POI05431 C COMPUTE THE COEFFICIENTS AN,BN,CN CORRESPONDING TO THE T DIRECTIONPOI05432 C POI05433 HDT = DELTAT/2. POI05434 TDT = DELTAT+DELTAT POI05435 DO 104 J=1,N POI05436 TEMP1 = 1./(T(J)*TDT) POI05437 TEMP2 = 1./((T(J)-HDT)*TDT) POI05438 TEMP3 = 1./((T(J)+HDT)*TDT) POI05439 AN(J) = TEMP1*TEMP2 POI05440 CN(J) = TEMP1*TEMP3 POI05441 BN(J) = -(AN(J)+CN(J)) POI05442 104 CONTINUE POI05443 C POI05444 C COMPUTE RIGHT SIDE OF EQUATION POI05445 C POI05446 DO 106 J=1,N POI05447 DO 105 I=1,M POI05448 Y(I,J) = 3.75*S(I)*T(J)*(S(I)**4+T(J)**4) POI05449 105 CONTINUE POI05450 106 CONTINUE POI05451 C POI05452 C INCLUDE NONHOMOGENEOUS BOUNDARY INTO RIGHT SIDE. NOTE THAT THE POI05453 C CORNER AT J=N,I=M INCLUDES CONTRIBUTIONS FROM BOTH BOUNDARIES. POI05454 C POI05455 DO 107 J=1,N POI05456 Y(M,J) = Y(M,J)-CM(M)*T(J)**5 POI05457 107 CONTINUE POI05458 DO 108 I=1,M POI05459 Y(I,N) = Y(I,N)-CN(N)*S(I)**5 POI05460 108 CONTINUE POI05461 C POI05462 109 CALL BLKTRI (IFLG,NP,N,AN,BN,CN,MP,M,AM,BM,CM,IDIMY,Y,IERROR,W) POI05463 IFLG = IFLG+1 POI05464 IF (IFLG-1) 109,109,110 POI05465 C POI05466 C COMPUTE DISCRETIZATION ERROR POI05467 C POI05468 110 ERR = 0. POI05469 DO 112 J=1,N POI05470 DO 111 I=1,M POI05471 Z = ABS(Y(I,J)-(S(I)*T(J))**5) POI05472 IF (Z .GT. ERR) ERR = Z POI05473 111 CONTINUE POI05474 112 CONTINUE POI05475 PRINT 1001 , IERROR,ERR POI05476 STOP POI05477 C POI05478 C POI05479 C POI05480 1001 FORMAT (///9H IERROR =,I3,10X,22HDISCRETIZATION ERROR =,E12.4) POI05481 C POI05482 END POI05483