SLATEC Routines --- HWSCSP ---
*DECK HWSCSP
SUBROUTINE HWSCSP (INTL, TS, TF, M, MBDCND, BDTS, BDTF, RS, RF, N,
+ NBDCND, BDRS, BDRF, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
C***BEGIN PROLOGUE HWSCSP
C***PURPOSE Solve a finite difference approximation to the modified
C Helmholtz equation in spherical coordinates assuming
C axisymmetry (no dependence on longitude).
C***LIBRARY SLATEC (FISHPACK)
C***CATEGORY I2B1A1A
C***TYPE SINGLE PRECISION (HWSCSP-S)
C***KEYWORDS ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, SPHERICAL
C***AUTHOR Adams, J., (NCAR)
C Swarztrauber, P. N., (NCAR)
C Sweet, R., (NCAR)
C***DESCRIPTION
C
C Subroutine HWSCSP 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
C This two dimensional modified Helmholtz equation results from
C the Fourier transform of the three dimensional Poisson equation
C
C * * * * * * * * * * On Input * * * * * * * * * *
C
C INTL
C = 0 On initial entry to HWSCSP 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 HWSCSP.
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
C * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
C
C If TF is equal to PI then it must be computed using the statement
C TF = PIMACH(DUM). This insures that TF in the users program is
C equal to PI in this program which permits several tests of the
C input parameters that otherwise would not be possible.
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(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.
C N must be greater than 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, HWSCSP 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 HWSCSP. 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. Its length can be computed from the formula below
C which depends on the value 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 L=2**(K+1) then W must be
C dimensioned at least (K-2)*L+K+5*(M+N)+MAX(2*N,6*M)+23
C
C **IMPORTANT** For purposes of checking, the required length
C of W is computed by HWSCSP and stored in W(1)
C in floating point format.
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. HWSCSP 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 = 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 7
C
C Since this is the only means of indicating a possibly incorrect
C call to HWSCSP, the user should test IERROR after a call.
C
C W
C Contains intermediate values that must not be destroyed if
C HWSCSP will be called again with INTL = 1. W(1) contains the
C number of locations which W must have.
C
C *Long Description:
C
C * * * * * * * Program Specifications * * * * * * * * * * * *
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 June 1979
C Revision
C
C Subprograms HWSCSP,HWSCS1,BLKTRI,BLKTR1,PROD,PRODP,CPROD,CPRODP
C Required ,COMBP,PPADD,PSGF,BSRH,PPSGF,PPSPF,TEVLS,INDXA,
C ,INDXB,INDXC,R1MACH
C
C Special
C Conditions
C
C Common CBLKT
C Blocks
C
C I/O NONE
C
C Precision Single
C
C Specialist Paul N Swarztrauber
C
C Language FORTRAN
C
C History Version 1 September 1973
C Version 2 April 1976
C Version 3 June 1979
C
C Algorithm The routine defines the finite difference
C equations, incorporates boundary data, and adjusts
C the right side of singular systems and then calls
C BLKTRI to solve the system.
C
C Space
C Required
C
C Portability American National Standards Institute FORTRAN.
C The machine accuracy is set using function R1MACH.
C
C Required NONE
C Resident
C Routines
C
C Reference Swarztrauber,P. and R. Sweet, 'Efficient FORTRAN
C Subprograms for The Solution Of Elliptic Equations'
C NCAR TN/IA-109, July, 1975, 138 pp.
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES P. N. Swarztrauber and R. Sweet, Efficient Fortran
C subprograms for the solution of elliptic equations,
C NCAR TN/IA-109, July 1975, 138 pp.
C***ROUTINES CALLED HWSCS1, PIMACH
C***REVISION HISTORY (YYMMDD)
C 801001 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE HWSCSP