SUBROUTINE HSTCSP (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N,NBDCND,BDC, 1 BDD,ELMBDA,F,IDIMF,PERTRB,IERROR,W) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F I S H P A C K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 4.0 , JUNE 1999) * C * * C * BY * C * * C * JOHN ADAMS, PAUL SWARZTRAUBER AND ROLAND SWEET * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C C DIMENSION OF BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N), C ARGUMENTS W(SEE ARGUMENT LIST) C C LATEST REVISION NOVEMBER 1988 C C PURPOSE SOLVES THE STANDARD FIVE-POINT FINITE C DIFFERENCE APPROXIMATION ON A STAGGERED C GRID TO THE MODIFIED HELMHOLTZ EQUATION IN C SPHERICAL COORDINATES ASSUMING AXISYMMETRY C (NO DEPENDENCE ON LONGITUDE). C C THE EQUATION IS C C (1/R**2)(D/DR)(R**2(DU/DR)) + C 1/(R**2*SIN(THETA))(D/DTHETA) C (SIN(THETA)(DU/DTHETA)) + C (LAMBDA/(R*SIN(THETA))**2)U = F(THETA,R) C C WHERE THETA IS COLATITUDE AND R IS THE C RADIAL COORDINATE. THIS TWO-DIMENSIONAL C MODIFIED HELMHOLTZ EQUATION RESULTS FROM C THE FOURIER TRANSFORM OF THE THREE- C DIMENSIONAL POISSON EQUATION. C C C USAGE CALL HSTCSP (INTL,A,B,M,MBDCND,BDA,BDB,C,D,N, C NBDCND,BDC,BDD,ELMBDA,F,IDIMF, C PERTRB,IERROR,W) C C ARGUMENTS C ON INPUT INTL C C = 0 ON INITIAL ENTRY TO HSTCSP OR IF ANY C OF THE ARGUMENTS C, D, N, OR NBDCND C ARE CHANGED FROM A PREVIOUS CALL C C = 1 IF C, D, N, AND NBDCND ARE ALL C UNCHANGED FROM PREVIOUS CALL TO HSTCSP C C NOTE: C A CALL WITH INTL = 0 TAKES APPROXIMATELY C 1.5 TIMES AS MUCH TIME AS A CALL WITH C INTL = 1. ONCE A CALL WITH INTL = 0 C HAS BEEN MADE THEN SUBSEQUENT SOLUTIONS C CORRESPONDING TO DIFFERENT F, BDA, BDB, C BDC, AND BDD CAN BE OBTAINED FASTER WITH C INTL = 1 SINCE INITIALIZATION IS NOT C REPEATED. C C A,B C THE RANGE OF THETA (COLATITUDE), C I.E. A .LE. THETA .LE. B. A C MUST BE LESS THAN B AND A MUST BE C NON-NEGATIVE. A AND B ARE IN RADIANS. C A = 0 CORRESPONDS TO THE NORTH POLE AND C B = PI CORRESPONDS TO THE SOUTH POLE. C C * * * IMPORTANT * * * C C IF B IS EQUAL TO PI, THEN B MUST BE C COMPUTED USING THE STATEMENT C B = PIMACH(DUM) C THIS INSURES THAT B IN THE USER'S PROGRAM C IS EQUAL TO PI IN THIS PROGRAM, PERMITTING C SEVERAL TESTS OF THE INPUT PARAMETERS THAT C OTHERWISE WOULD NOT BE POSSIBLE. C C * * * * * * * * * * * * C C M C THE NUMBER OF GRID POINTS IN THE INTERVAL C (A,B). THE GRID POINTS IN THE THETA- C DIRECTION ARE GIVEN BY C THETA(I) = A + (I-0.5)DTHETA C FOR I=1,2,...,M WHERE DTHETA =(B-A)/M. C M MUST BE GREATER THAN 4. C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS C AT THETA = A AND THETA = B. C C = 1 IF THE SOLUTION IS SPECIFIED AT C THETA = A AND THETA = B. C (SEE NOTES 1, 2 BELOW) C C = 2 IF THE SOLUTION IS SPECIFIED AT C THETA = A AND THE DERIVATIVE OF THE C SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = B C (SEE NOTES 1, 2 BELOW). C C = 3 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED C AT THETA = A (SEE NOTES 1, 2 BELOW) C AND THETA = B. C C = 4 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED AT C THETA = A (SEE NOTES 1, 2 BELOW) AND C THE SOLUTION IS SPECIFIED AT THETA = B. C C = 5 IF THE SOLUTION IS UNSPECIFIED AT C THETA = A = 0 AND THE SOLUTION IS C SPECIFIED AT THETA = B. C (SEE NOTE 2 BELOW) C C = 6 IF THE SOLUTION IS UNSPECIFIED AT C THETA = A = 0 AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = B C (SEE NOTE 2 BELOW). C C = 7 IF THE SOLUTION IS SPECIFIED AT C THETA = A AND THE SOLUTION IS C UNSPECIFIED AT THETA = B = PI. C C = 8 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED AT C THETA = A (SEE NOTE 1 BELOW) C AND THE SOLUTION IS UNSPECIFIED AT C THETA = B = PI. C C = 9 IF THE SOLUTION IS UNSPECIFIED AT C THETA = A = 0 AND THETA = B = PI. C C NOTE 1: C IF A = 0, DO NOT USE MBDCND = 1,2,3,4,7 C OR 8, BUT INSTEAD USE MBDCND = 5, 6, OR 9. C C NOTE 2: C IF B = PI, DO NOT USE MBDCND = 1,2,3,4,5, C OR 6, BUT INSTEAD USE MBDCND = 7, 8, OR 9. C C NOTE 3: C WHEN A = 0 AND/OR B = PI THE ONLY C MEANINGFUL BOUNDARY CONDITION IS C DU/DTHETA = 0. SEE D. GREENSPAN, C 'NUMERICAL ANALYSIS OF ELLIPTIC C BOUNDARY VALUE PROBLEMS,' C HARPER AND ROW, 1965, CHAPTER 5.) C C BDA C A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT C SPECIFIES THE BOUNDARY VALUES (IF ANY) OF C THE SOLUTION AT THETA = A. C C WHEN MBDCND = 1, 2, OR 7, C BDA(J) = U(A,R(J)), J=1,2,...,N. C C WHEN MBDCND = 3, 4, OR 8, C BDA(J) = (D/DTHETA)U(A,R(J)), J=1,2,...,N. C C WHEN MBDCND HAS ANY OTHER VALUE, BDA IS A C DUMMY VARIABLE. C C BDB C A ONE-DIMENSIONAL ARRAY OF LENGTH N THAT C SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT THETA = B. C C WHEN MBDCND = 1, 4, OR 5, C BDB(J) = U(B,R(J)), J=1,2,...,N. C C WHEN MBDCND = 2,3, OR 6, C BDB(J) = (D/DTHETA)U(B,R(J)), J=1,2,...,N. C C WHEN MBDCND HAS ANY OTHER VALUE, BDB IS C A DUMMY VARIABLE. C C C,D C THE RANGE OF R , I.E. C .LE. R .LE. D. C C MUST BE LESS THAN D AND NON-NEGATIVE. C C N C THE NUMBER OF UNKNOWNS IN THE INTERVAL C (C,D). THE UNKNOWNS IN THE R-DIRECTION C ARE GIVEN BY R(J) = C + (J-0.5)DR, C J=1,2,...,N, WHERE DR = (D-C)/N. C N MUST BE GREATER THAN 4. C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITIONS C AT R = C AND R = D. C C C = 1 IF THE SOLUTION IS SPECIFIED AT C R = C AND R = D. C C = 2 IF THE SOLUTION IS SPECIFIED AT C R = C AND THE DERIVATIVE OF THE C SOLUTION WITH RESPECT TO R IS C SPECIFIED AT R = D. (SEE NOTE 1 BELOW) C C = 3 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO R IS SPECIFIED AT C R = C AND R = D. C C = 4 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO R IS C SPECIFIED AT R = C AND THE SOLUTION C IS SPECIFIED AT R = D. C C = 5 IF THE SOLUTION IS UNSPECIFIED AT C R = C = 0 (SEE NOTE 2 BELOW) AND THE C SOLUTION IS SPECIFIED AT R = D. C C = 6 IF THE SOLUTION IS UNSPECIFIED AT C R = C = 0 (SEE NOTE 2 BELOW) C AND THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO R IS SPECIFIED AT C R = D. C C NOTE 1: C IF C = 0 AND MBDCND = 3,6,8 OR 9, THE C SYSTEM OF EQUATIONS TO BE SOLVED IS C SINGULAR. THE UNIQUE SOLUTION IS C DETERMINED BY EXTRAPOLATION TO THE C SPECIFICATION OF U(THETA(1),C). C BUT IN THESE CASES THE RIGHT SIDE OF THE C SYSTEM WILL BE PERTURBED BY THE CONSTANT C PERTRB. C C NOTE 2: C NBDCND = 5 OR 6 CANNOT BE USED WITH C MBDCND =1, 2, 4, 5, OR 7 C (THE FORMER INDICATES THAT THE SOLUTION IS C UNSPECIFIED AT R = 0; THE LATTER INDICATES C SOLUTION IS SPECIFIED). C USE INSTEAD NBDCND = 1 OR 2. C C BDC C A ONE DIMENSIONAL ARRAY OF LENGTH M THAT C SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT R = C. WHEN NBDCND = 1 OR 2, C BDC(I) = U(THETA(I),C), I=1,2,...,M. C C WHEN NBDCND = 3 OR 4, C BDC(I) = (D/DR)U(THETA(I),C), I=1,2,...,M. C C WHEN NBDCND HAS ANY OTHER VALUE, BDC IS C A DUMMY VARIABLE. C C BDD C A ONE-DIMENSIONAL ARRAY OF LENGTH M THAT C SPECIFIES THE BOUNDARY VALUES OF THE C SOLUTION AT R = D. WHEN NBDCND = 1 OR 4, C BDD(I) = U(THETA(I),D) , I=1,2,...,M. C C WHEN NBDCND = 2 OR 3, C BDD(I) = (D/DR)U(THETA(I),D), I=1,2,...,M. C C WHEN NBDCND HAS ANY OTHER VALUE, BDD IS C A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE MODIFIED C HELMHOLTZ EQUATION. IF LAMBDA IS GREATER C THAN 0, A SOLUTION MAY NOT EXIST. C HOWEVER, HSTCSP WILL ATTEMPT TO FIND A C SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE C VALUES OF THE RIGHT SIDE OF THE MODIFIED C HELMHOLTZ EQUATION. FOR I=1,2,...,M AND C J=1,2,...,N C C F(I,J) = F(THETA(I),R(J)) . C C F MUST BE DIMENSIONED AT LEAST M X N. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY C F AS IT APPEARS IN THE PROGRAM CALLING C HSTCSP. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. C IDIMF MUST BE AT LEAST M. C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE C PROVIDED BY THE USER FOR WORK SPACE. C WITH K = INT(LOG2(N))+1 AND L = 2**(K+1), C W MAY REQUIRE UP TO C (K-2)*L+K+MAX(2N,6M)+4(N+M)+5 LOCATIONS. C THE ACTUAL NUMBER OF LOCATIONS USED IS C COMPUTED BY HSTCSP AND IS RETURNED IN THE C LOCATION W(1). C C C ON OUTPUT F C CONTAINS THE SOLUTION U(I,J) OF THE FINITE C DIFFERENCE APPROXIMATION FOR THE GRID POINT C (THETA(I),R(J)) FOR I=1,2,..,M, J=1,2,...,N. C C PERTRB C IF A COMBINATION OF PERIODIC, DERIVATIVE, C OR UNSPECIFIED BOUNDARY CONDITIONS IS C SPECIFIED FOR A POISSON EQUATION C (LAMBDA = 0), A SOLUTION MAY NOT EXIST. C PERTRB IS A CONSTANT, CALCULATED AND C SUBTRACTED FROM F, WHICH ENSURES THAT A C SOLUTION EXISTS. HSTCSP THEN COMPUTES THIS C SOLUTION, WHICH IS A LEAST SQUARES SOLUTION C TO THE ORIGINAL APPROXIMATION. C THIS SOLUTION PLUS ANY CONSTANT IS ALSO C A SOLUTION; HENCE, THE SOLUTION IS NOT C UNIQUE. THE VALUE OF PERTRB SHOULD BE C SMALL COMPARED TO THE RIGHT SIDE F. C OTHERWISE, A SOLUTION IS OBTAINED TO AN C ESSENTIALLY DIFFERENT PROBLEM. C THIS COMPARISON SHOULD ALWAYS BE MADE TO C INSURE THAT A MEANINGFUL SOLUTION HAS BEEN C OBTAINED. C C IERROR C AN ERROR FLAG THAT INDICATES INVALID INPUT C PARAMETERS. EXCEPT FOR NUMBERS 0 AND 10, C A SOLUTION IS NOT ATTEMPTED. C C = 0 NO ERROR C C = 1 A .LT. 0 OR B .GT. PI C C = 2 A .GE. B C C = 3 MBDCND .LT. 1 OR MBDCND .GT. 9 C C = 4 C .LT. 0 C C = 5 C .GE. D C C = 6 NBDCND .LT. 1 OR NBDCND .GT. 6 C C = 7 N .LT. 5 C C = 8 NBDCND = 5 OR 6 AND C MBDCND = 1, 2, 4, 5, OR 7 C C = 9 C .GT. 0 AND NBDCND .GE. 5 C C = 10 ELMBDA .GT. 0 C C = 11 IDIMF .LT. M C C = 12 M .LT. 5 C C = 13 A = 0 AND MBDCND =1,2,3,4,7 OR 8 C C = 14 B = PI AND MBDCND .LE. 6 C C = 15 A .GT. 0 AND MBDCND = 5, 6, OR 9 C C = 16 B .LT. PI AND MBDCND .GE. 7 C C = 17 LAMBDA .NE. 0 AND NBDCND .GE. 5 C C SINCE THIS IS THE ONLY MEANS OF INDICATING C A POSSIBLY INCORRECT CALL TO HSTCSP, C THE USER SHOULD TEST IERROR AFTER THE CALL. C C W C W(1) CONTAINS THE REQUIRED LENGTH OF W. C ALSO W CONTAINS INTERMEDIATE VALUES THAT C MUST NOT BE DESTROYED IF HSTCSP WILL BE C CALLED AGAIN WITH INTL = 1. C C I/O NONE C C PRECISION SINGLE C C REQUIRED LIBRARY BLKTRI AND COMF FROM FISHPAK C FILES C C LANGUAGE FORTRAN C C HISTORY WRITTEN BY ROLAND SWEET AT NCAR IN 1977. C RELEASED ON NCAR'S PUBLIC SOFTWARE LIBRARIES C IN JANUARY 1980. C C PORTABILITY FORTRAN 77 C C ALGORITHM THIS SUBROUTINE DEFINES THE FINITE-DIFFERENCE C EQUATIONS, INCORPORATES BOUNDARY DATA, ADJUSTS C THE RIGHT SIDE WHEN THE SYSTEM IS SINGULAR C AND CALLS BLKTRI WHICH SOLVES THE LINEAR C SYSTEM OF EQUATIONS. C C C TIMING FOR LARGE M AND N, THE OPERATION COUNT IS C ROUGHLY PROPORTIONAL TO M*N*LOG2(N). THE C TIMING ALSO DEPENDS ON INPUT PARAMETER INTL. C C ACCURACY THE SOLUTION PROCESS EMPLOYED RESULTS IN C A LOSS OF NO MORE THAN FOUR SIGNIFICANT C DIGITS FOR N AND M AS LARGE AS 64. C MORE DETAILED INFORMATION ABOUT ACCURACY C CAN BE FOUND IN THE DOCUMENTATION FOR C SUBROUTINE BLKTRI WHICH IS THE ROUTINE C SOLVES THE FINITE DIFFERENCE EQUATIONS. C C REFERENCES P.N. SWARZTRAUBER, "A DIRECT METHOD FOR C THE DISCRETE SOLUTION OF SEPARABLE ELLIPTIC C EQUATIONS", C SIAM J. NUMER. ANAL. 11(1974), PP. 1136-1150. C C U. SCHUMANN AND R. SWEET, "A DIRECT METHOD FOR C THE SOLUTION OF POISSON'S EQUATION WITH NEUMANN C BOUNDARY CONDITIONS ON A STAGGERED GRID OF C ARBITRARY SIZE," J. COMP. PHYS. 20(1976), C PP. 171-182. C***********************************************************************