SUBROUTINE HWSSSP (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF,N,NBDCND,BDPS, 1 BDPF,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 BDTS(N+1), BDTF(N+1), BDPS(M+1), BDPF(M+1), C ARGUMENTS F(IDIMF,N+1), W(SEE ARGUMENT LIST) C C LATEST REVISION NOVEMBER 1988 C C PURPOSE SOLVES A FINITE DIFFERENCE APPROXIMATION TO C THE HELMHOLTZ EQUATION IN SPHERICAL C COORDINATES AND ON THE SURFACE OF THE UNIT C SPHERE (RADIUS OF 1). THE EQUATION IS C C (1/SIN(THETA))(D/DTHETA)(SIN(THETA) C (DU/DTHETA)) + (1/SIN(THETA)**2)(D/DPHI) C (DU/DPHI) + LAMBDA*U = F(THETA,PHI) C C WHERE THETA IS COLATITUDE AND PHI IS C LONGITUDE. C C USAGE CALL HWSSSP (TS,TF,M,MBDCND,BDTS,BDTF,PS,PF, C N,NBDCND,BDPS,BDPF,ELMBDA,F, C IDIMF,PERTRB,IERROR,W) C C ARGUMENTS C ON INPUT TS,TF C C THE RANGE OF THETA (COLATITUDE), I.E., C TS .LE. THETA .LE. TF. TS MUST BE LESS C THAN TF. TS AND TF ARE IN RADIANS. C A TS OF ZERO CORRESPONDS TO THE NORTH C POLE AND A TF OF PI CORRESPONDS TO C THE SOUTH POLE. C C * * * IMPORTANT * * * C C IF TF IS EQUAL TO PI THEN IT MUST BE C COMPUTED USING THE STATEMENT C TF = PIMACH(DUM). THIS INSURES THAT TF C IN THE USER'S PROGRAM IS EQUAL TO PI IN C THIS PROGRAM WHICH PERMITS SEVERAL TESTS C OF THE INPUT PARAMETERS THAT OTHERWISE C WOULD NOT BE POSSIBLE. C C C M C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (TS,TF) IS SUBDIVIDED. C HENCE, THERE WILL BE M+1 GRID POINTS IN THE C THETA-DIRECTION GIVEN BY C THETA(I) = (I-1)DTHETA+TS FOR C I = 1,2,...,M+1, WHERE C DTHETA = (TF-TS)/M IS THE PANEL WIDTH. C M MUST BE GREATER THAN 5 C C MBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION C AT THETA = TS AND THETA = TF. C C = 1 IF THE SOLUTION IS SPECIFIED AT C THETA = TS AND THETA = TF. C = 2 IF THE SOLUTION IS SPECIFIED AT C THETA = TS AND THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO THETA IS C SPECIFIED AT THETA = TF C (SEE NOTE 2 BELOW). C = 3 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED C SPECIFIED AT THETA = TS AND C THETA = TF (SEE NOTES 1,2 BELOW). C = 4 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED C AT THETA = TS (SEE NOTE 1 BELOW) C AND THE SOLUTION IS SPECIFIED AT C THETA = TF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT C THETA = TS = 0 AND THE SOLUTION C IS SPECIFIED AT THETA = TF. C = 6 IF THE SOLUTION IS UNSPECIFIED AT C THETA = TS = 0 AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO THETA C IS SPECIFIED AT THETA = TF C (SEE NOTE 2 BELOW). C = 7 IF THE SOLUTION IS SPECIFIED AT C THETA = TS AND THE SOLUTION IS C IS UNSPECIFIED AT THETA = TF = PI. C = 8 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO THETA IS SPECIFIED C AT THETA = TS (SEE NOTE 1 BELOW) AND C THE SOLUTION IS UNSPECIFIED AT C THETA = TF = PI. C = 9 IF THE SOLUTION IS UNSPECIFIED AT C THETA = TS = 0 AND THETA = TF = PI. C C NOTES: C IF TS = 0, DO NOT USE MBDCND = 3,4, OR 8, C BUT INSTEAD USE MBDCND = 5,6, OR 9 . C C IF TF = PI, DO NOT USE MBDCND = 2,3, OR 6, C BUT INSTEAD USE MBDCND = 7,8, OR 9 . C C BDTS C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 THAT C SPECIFIES THE VALUES OF THE DERIVATIVE OF C 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)), C J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTS IS C A DUMMY VARIABLE. C C BDTF C A ONE-DIMENSIONAL ARRAY OF LENGTH N+1 C THAT SPECIFIES THE VALUES OF THE DERIVATIVE C 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)), C J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTF IS C A DUMMY VARIABLE. C C PS,PF C THE RANGE OF PHI (LONGITUDE), I.E., C PS .LE. PHI .LE. PF. PS MUST BE LESS C THAN PF. PS AND PF ARE IN RADIANS. C IF PS = 0 AND PF = 2*PI, PERIODIC C BOUNDARY CONDITIONS ARE USUALLY PRESCRIBED. C C * * * IMPORTANT * * * C C IF PF IS EQUAL TO 2*PI THEN IT MUST BE C COMPUTED USING THE STATEMENT C PF = 2.*PIMACH(DUM). THIS INSURES THAT C PF IN THE USERS PROGRAM IS EQUAL TO C 2*PI IN THIS PROGRAM WHICH PERMITS TESTS C OF THE INPUT PARAMETERS THAT OTHERWISE C WOULD NOT BE POSSIBLE. C C N C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (PS,PF) IS SUBDIVIDED. C HENCE, THERE WILL BE N+1 GRID POINTS C IN THE PHI-DIRECTION GIVEN BY C PHI(J) = (J-1)DPHI+PS FOR C J = 1,2,...,N+1, WHERE C DPHI = (PF-PS)/N IS THE PANEL WIDTH. C N MUST BE GREATER THAN 4 C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION C AT PHI = PS AND PHI = PF. C C = 0 IF THE SOLUTION IS PERIODIC IN PHI, C I.U., U(I,J) = U(I,N+J). C = 1 IF THE SOLUTION IS SPECIFIED AT C PHI = PS AND PHI = PF C (SEE NOTE BELOW). C = 2 IF THE SOLUTION IS SPECIFIED AT C PHI = PS (SEE NOTE BELOW) C AND THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO PHI IS SPECIFIED C AT PHI = PF. C = 3 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO PHI IS SPECIFIED C AT PHI = PS AND PHI = PF. C = 4 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO PHI IS SPECIFIED C AT PS AND THE SOLUTION IS SPECIFIED C AT PHI = PF C C NOTE: C NBDCND = 1,2, OR 4 CANNOT BE USED WITH C MBDCND = 5,6,7,8, OR 9. THE FORMER INDICATES C THAT THE SOLUTION IS SPECIFIED AT A POLE, THE C LATTER INDICATES THAT THE SOLUTION IS NOT C SPECIFIED. USE INSTEAD MBDCND = 1 OR 2. C C BDPS C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT C SPECIFIES THE VALUES OF THE DERIVATIVE C 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), C I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDPS IS C A DUMMY VARIABLE. C C BDPF C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT C SPECIFIES THE VALUES OF THE DERIVATIVE C 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), C I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDPF IS C A DUMMY VARIABLE. C C ELMBDA C THE CONSTANT LAMBDA IN THE HELMHOLTZ C EQUATION. IF LAMBDA .GT. 0, A SOLUTION C MAY NOT EXIST. HOWEVER, HWSSSP WILL C ATTEMPT TO FIND A SOLUTION. C C F C A TWO-DIMENSIONAL ARRAY THAT SPECIFIES THE C VALUE OF THE RIGHT SIDE OF THE HELMHOLTZ C EQUATION AND BOUNDARY VALUES (IF ANY). C F MUST BE DIMENSIONED AT LEAST (M+1)*(N+1). C C ON THE INTERIOR, F IS DEFINED AS FOLLOWS: C FOR I = 2,3,...,M AND J = 2,3,...,N C F(I,J) = F(THETA(I),PHI(J)). C C ON THE BOUNDARIES F IS DEFINED AS FOLLOWS: C FOR J = 1,2,...,N+1 AND I = 1,2,...,M+1 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)) 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) C 3 F(THETA(I),PS) F(THETA(I),PF) C 4 F(THETA(I),PS) U(THETA(I),PF) C C NOTE: C IF THE TABLE CALLS FOR BOTH THE SOLUTION U C AND THE RIGHT SIDE F AT A CORNER THEN THE C SOLUTION MUST BE SPECIFIED. C C IDIMF C THE ROW (OR FIRST) DIMENSION OF THE ARRAY C F AS IT APPEARS IN THE PROGRAM CALLING C HWSSSP. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. IDIMF MUST BE C AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE C PROVIDED BY THE USER FOR WORK SPACE. C W MAY REQUIRE UP TO C 4*(N+1)+(16+INT(LOG2(N+1)))(M+1) LOCATIONS C THE ACTUAL NUMBER OF LOCATIONS USED IS C COMPUTED BY HWSSSP AND IS OUTPUT IN C LOCATION W(1). INT( ) DENOTES THE C FORTRAN INTEGER FUNCTION. 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),PHI(J)), I = 1,2,...,M+1 AND C J = 1,2,...,N+1 . C C PERTRB C IF ONE SPECIFIES A COMBINATION OF PERIODIC, C DERIVATIVE OR UNSPECIFIED BOUNDARY C CONDITIONS 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. HWSSSP THEN COMPUTES C THIS SOLUTION, WHICH IS A LEAST SQUARES C SOLUTION TO THE ORIGINAL APPROXIMATION. C THIS SOLUTION IS NOT UNIQUE AND IS C UNNORMALIZED. THE VALUE OF PERTRB SHOULD C BE SMALL COMPARED TO THE RIGHT SIDE F. C OTHERWISE , A SOLUTION IS OBTAINED TO AN C ESSENTIALLY DIFFERENT PROBLEM. THIS C COMPARISON SHOULD ALWAYS BE MADE TO INSURE C 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 8, C 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.LT.0 OR PS.GT.PI+PI C = 5 PS.GE.PF C = 6 N.LT.5 C = 7 M.LT.5 C = 8 NBDCND.LT.0 OR NBDCND.GT.4 C = 9 ELMBDA.GT.0 C = 10 IDIMF.LT.M+1 C = 11 NBDCND EQUALS 1,2 OR 4 AND MBDCND.GE.5 C = 12 TS.EQ.0 AND MBDCND EQUALS 3,4 OR 8 C = 13 TF.EQ.PI AND MBDCND EQUALS 2,3 OR 6 C = 14 MBDCND EQUALS 5,6 OR 9 AND TS.NE.0 C = 15 MBDCND.GE.7 AND TF.NE.PI C C SINCE THIS IS THE ONLY MEANS OF INDICATING C A POSSIBLY INCORRECT CALL TO HWSSSP, THE C USER SHOULD TEST IERROR AFTER A CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT C BE DESTROYED IF HWSSSP WILL BE CALLED AGAIN C WITH INTL = 1. W(1) CONTAINS THE REQUIRED C LENGTH OF W . C C SPECIAL CONDITIONS NONE C C I/O NONE C C PRECISION SINGLE C C REQUIRED LIBRARY GENBUN, GNBNAUX, AND COMF C FILES FROM FISHPAK C C LANGUAGE FORTRAN C C HISTORY WRITTEN BY ROLAND SWEET AT NCAR IN THE LATE C 1970'S. RELEASED ON NCAR'S PUBLIC SOFTWARE C LIBRARIES IN JANUARY 1980. C C PORTABILITY FORTRAN 77 C C ALGORITHM THE ROUTINE DEFINES THE FINITE DIFFERENCE C EQUATIONS, INCORPORATES BOUNDARY DATA, AND C ADJUSTS THE RIGHT SIDE OF SINGULAR SYSTEMS C AND THEN CALLS GENBUN TO SOLVE THE SYSTEM. C C TIMING FOR LARGE M AND N, THE OPERATION COUNT C IS ROUGHLY PROPORTIONAL TO C M*N*(LOG2(N) C BUT ALSO DEPENDS ON INPUT PARAMETERS NBDCND C AND MBDCND. C C ACCURACY THE SOLUTION PROCESS EMPLOYED RESULTS IN A LOSS C OF NO MORE THAN THREE SIGNIFICANT DIGITS FOR N C AND M AS LARGE AS 64. MORE DETAILS ABOUT C ACCURACY CAN BE FOUND IN THE DOCUMENTATION FOR C SUBROUTINE GENBUN WHICH IS THE ROUTINE THAT C SOLVES THE FINITE DIFFERENCE EQUATIONS. C C REFERENCES P. N. SWARZTRAUBER, "THE DIRECT SOLUTION OF C THE DISCRETE POISSON EQUATION ON THE SURFACE OF C A SPHERE", S.I.A.M. J. NUMER. ANAL.,15(1974), C PP 212-215. C C SWARZTRAUBER,P. AND R. SWEET, "EFFICIENT C FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC EQUATIONS", NCAR TN/IA-109, JULY, C 1975, 138 PP. C***********************************************************************