SUBROUTINE HWSCSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF,RS,RF,N,NBDCND, 1 BDRS,BDRF,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), BDRS(M+1), BDRF(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 C TO THE MODIFIED HELMHOLTZ EQUATION IN C SPHERICAL COORDINATES ASSUMING AXISYMMETRY C (NO DEPENDENCE ON LONGITUDE). THE EQUATION C IS C C (1/R**2)(D/DR)((R**2)(D/DR)U) + C C (1/(R**2)SIN(THETA))(D/DTHETA) C C (SIN(THETA)(D/DTHETA)U) + C C (LAMBDA/(RSIN(THETA))**2)U = F(THETA,R). C C THIS TWO DIMENSIONAL MODIFIED HELMHOLTZ C EQUATION RESULTS FROM THE FOURIER TRANSFORM C OF THE THREE DIMENSIONAL POISSON EQUATION. C C USAGE CALL HWSCSP (INTL,TS,TF,M,MBDCND,BDTS,BDTF, C RS,RF,N,NBDCND,BDRS,BDRF,ELMBDA, C F,IDIMF,PERTRB,IERROR,W) C C ARGUMENTS C ON INPUT INTL C = 0 ON INITIAL ENTRY TO HWSCSP OR IF ANY C OF THE ARGUMENTS RS, RF, N, NBDCND C ARE CHANGED FROM A PREVIOUS CALL. C = 1 IF RS, RF, N, NBDCND ARE ALL UNCHANGED C FROM PREVIOUS CALL TO HWSCSP. 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, BDTS, BDTF, C BDRS, BDRF CAN BE OBTAINED FASTER WITH C INTL = 1 SINCE INITIALIZATION IS NOT C REPEATED. C C TS,TF 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. A TS OF C ZERO CORRESPONDS TO THE NORTH POLE AND A C TF OF PI CORRESPONDS TO 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 M C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (TS,TF) IS SUBDIVIDED. C HENCE, THERE WILL BE M+1 GRID POINTS C IN THE THETA-DIRECTION GIVEN BY C THETA(K) = (I-1)DTHETA+TS FOR C I = 1,2,...,M+1, WHERE DTHETA = (TF-TS)/M C IS THE PANEL WIDTH. 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 THE C 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 AT THETA = TS AND THETA = TF C (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) AND C SOLUTION IS SPECIFIED AT THETA = TF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT C THETA = TS = 0 AND THE SOLUTION IS C 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 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) C AND 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 NOTE 1: C IF TS = 0, DO NOT USE MBDCND = 3,4, OR 8, C BUT INSTEAD USE MBDCND = 5,6, OR 9 . C C NOTE 2: 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,R(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 THAT C SPECIFIES THE VALUES OF THE DERIVATIVE OF C 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)), C J = 1,2,...,N+1 . C C WHEN MBDCND HAS ANY OTHER VALUE, BDTF IS C A DUMMY VARIABLE. C C RS,RF C THE RANGE OF R, I.E., RS .LE. R .LT. RF. C RS MUST BE LESS THAN RF. RS MUST BE C NON-NEGATIVE. C C N C THE NUMBER OF PANELS INTO WHICH THE C INTERVAL (RS,RF) IS SUBDIVIDED. C HENCE, THERE WILL BE N+1 GRID POINTS IN THE C R-DIRECTION GIVEN BY R(J) = (J-1)DR+RS C FOR J = 1,2,...,N+1, WHERE DR = (RF-RS)/N C IS THE PANEL WIDTH. C N MUST BE GREATER THAN 2 C C NBDCND C INDICATES THE TYPE OF BOUNDARY CONDITION C AT R = RS AND R = RF. C C = 1 IF THE SOLUTION IS SPECIFIED AT C R = RS AND R = RF. C = 2 IF THE SOLUTION IS SPECIFIED AT C R = RS AND THE DERIVATIVE C OF THE SOLUTION WITH RESPECT TO R C IS SPECIFIED AT R = RF. C = 3 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO R IS SPECIFIED AT C R = RS AND R = RF. C = 4 IF THE DERIVATIVE OF THE SOLUTION C WITH RESPECT TO R IS SPECIFIED AT C RS AND THE SOLUTION IS SPECIFIED AT C R = RF. C = 5 IF THE SOLUTION IS UNSPECIFIED AT C R = RS = 0 (SEE NOTE BELOW) AND THE C SOLUTION IS SPECIFIED AT R = RF. C = 6 IF THE SOLUTION IS UNSPECIFIED AT C R = RS = 0 (SEE NOTE BELOW) AND THE C DERIVATIVE OF THE SOLUTION WITH C RESPECT TO R IS SPECIFIED AT R = RF. C C NOTE: C NBDCND = 5 OR 6 CANNOT BE USED WITH C MBDCND = 1,2,4,5, OR 7. THE FORMER C INDICATES THAT THE SOLUTION IS UNSPECIFIED C AT R = 0, THE LATTER INDICATES THAT THE C SOLUTION IS SPECIFIED). C USE INSTEAD NBDCND = 1 OR 2 . C C BDRS C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 THAT C SPECIFIES THE VALUES OF THE DERIVATIVE OF C THE SOLUTION WITH RESPECT TO R AT R = RS. C C WHEN NBDCND = 3 OR 4, C BDRS(I) = (D/DR)U(THETA(I),RS), C I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDRS IS C A DUMMY VARIABLE. C C BDRF C A ONE-DIMENSIONAL ARRAY OF LENGTH M+1 C THAT SPECIFIES THE VALUES OF THE C DERIVATIVE OF THE SOLUTION WITH RESPECT C TO R AT R = RF. C C WHEN NBDCND = 2,3, OR 6, C BDRF(I) = (D/DR)U(THETA(I),RF), C I = 1,2,...,M+1 . C C WHEN NBDCND HAS ANY OTHER VALUE, BDRF 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, HWSCSP WILL C ATTEMPT TO FIND A SOLUTION. IF NBDCND = 5 C OR 6 OR MBDCND = 5,6,7,8, OR 9, ELMBDA C MUST BE ZERO. C C F C A TWO-DIMENSIONAL ARRAY, OF DIMENSION AT C LEAST (M+1)*(N+1), SPECIFYING VALUES OF THE C RIGHT SIDE OF THE HELMHOLTZ EQUATION AND C BOUNDARY VALUES (IF ANY). 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),R(J)). C C ON THE BOUNDARIES, F IS DEFINED AS FOLLOWS: C FOR J=1,2,...,N+1, I=1,2,...,M+1, 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)) 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) C 5 F(TS,0) U(THETA(I),RF) C 6 F(TS,0) F(THETA(I),RF) C C NOTE: C IF THE TABLE CALLS FOR BOTH THE SOLUTION C U AND THE RIGHT SIDE F AT A CORNER THEN C THE 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 HWSCSP. THIS PARAMETER IS USED TO SPECIFY C THE VARIABLE DIMENSION OF F. IDIMF MUST C BE AT LEAST M+1 . C C W C A ONE-DIMENSIONAL ARRAY THAT MUST BE C PROVIDED BY THE USER FOR WORK SPACE. C ITS LENGTH CAN BE COMPUTED FROM THE C FORMULA BELOW WHICH DEPENDS ON THE VALUE C OF NBDCND C C IF NBDCND=2,4 OR 6 DEFINE NUNK=N C IF NBDCND=1 OR 5 DEFINE NUNK=N-1 C IF NBDCND=3 DEFINE NUNK=N+1 C C NOW SET K=INT(LOG2(NUNK))+1 AND C L=2**(K+1) THEN W MUST BE DIMENSIONED C AT LEAST (K-2)*L+K+5*(M+N)+MAX(2*N,6*M)+23 C C **IMPORTANT** C FOR PURPOSES OF CHECKING, THE REQUIRED C LENGTH OF W IS COMPUTED BY HWSCSP AND C STORED IN W(1) IN FLOATING POINT FORMAT. 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)), I = 1,2,...,M+1, C J = 1,2,...,N+1 . C C PERTRB C IF A COMBINATION OF PERIODIC OR DERIVATIVE C BOUNDARY CONDITIONS IS SPECIFIED FOR A C POISSON EQUATION (LAMBDA = 0), A SOLUTION C MAY NOT EXIST. PERTRB IS A CONSTANT, C CALCULATED AND SUBTRACTED FROM F, WHICH C ENSURES THAT A SOLUTION EXISTS. HWSCSP C THEN COMPUTES THIS SOLUTION, WHICH IS A C LEAST SQUARES SOLUTION TO THE ORIGINAL C APPROXIMATION. THIS SOLUTION IS NOT UNIQUE C AND IS UNNORMALIZED. THE VALUE OF PERTRB C SHOULD BE SMALL COMPARED TO THE RIGHT SIDE C F. OTHERWISE , A SOLUTION IS OBTAINED TO C AN ESSENTIALLY DIFFERENT PROBLEM. THIS C COMPARISON SHOULD ALWAYS BE MADE TO INSURE C THAT A MEANINGFUL SOLUTION HAS BEEN 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 = 1 TS.LT.0. OR TF.GT.PI C = 2 TS.GE.TF C = 3 M.LT.5 C = 4 MBDCND.LT.1 OR MBDCND.GT.9 C = 5 RS.LT.0 C = 6 RS.GE.RF C = 7 N.LT.5 C = 8 NBDCND.LT.1 OR NBDCND.GT.6 C = 9 ELMBDA.GT.0 C = 10 IDIMF.LT.M+1 C = 11 ELMBDA.NE.0 AND MBDCND.GE.5 C = 12 ELMBDA.NE.0 AND NBDCND EQUALS 5 OR 6 C = 13 MBDCND EQUALS 5,6 OR 9 AND TS.NE.0 C = 14 MBDCND.GE.7 AND TF.NE.PI C = 15 TS.EQ.0 AND MBDCND EQUALS 3,4 OR 8 C = 16 TF.EQ.PI AND MBDCND EQUALS 2,3 OR 6 C = 17 NBDCND.GE.5 AND RS.NE.0 C = 18 NBDCND.GE.5 AND MBDCND EQUALS 1,2,4,5 OR C C SINCE THIS IS THE ONLY MEANS OF INDICATING C A POSSLIBY INCORRECT CALL TO HWSCSP, THE C USER SHOULD TEST IERROR AFTER A CALL. C C W C CONTAINS INTERMEDIATE VALUES THAT MUST NOT C BE DESTROYED IF HWSCSP WILL BE CALLED AGAIN C WITH INTL = 1. W(1) CONTAINS THE NUMBER C OF LOCATIONS WHICH W MUST HAVE C C SPECIAL CONDITIONS NONE 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 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 BLKTRI TO SOLVE THE SYSTEM. C C REFERENCES SWARZTRAUBER,P. AND R. SWEET, "EFFICIENT C FORTRAN SUBPROGRAMS FOR THE SOLUTION OF C ELLIPTIC EQUATIONS" C NCAR TN/IA-109, JULY, 1975, 138 PP. C***********************************************************************