# hwsssp.f

          SUBROUTINE HWSSSP (TS, TF, M, MBDCND, BDTS, BDTF, PS, PF, N,
         +   NBDCND, BDPS, BDPF, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
    C***BEGIN PROLOGUE  HWSSSP
    C***PURPOSE  Solve a finite difference approximation to the Helmholtz
    C            equation in spherical coordinates and on the surface of the
    C            unit sphere (radius of 1).
    C***LIBRARY   SLATEC (FISHPACK)
    C***CATEGORY  I2B1A1A
    C***TYPE      SINGLE PRECISION (HWSSSP-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 HWSSSP 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)(dU/dTHETA))
    C
    C             + (1/SIN(THETA)**2)(d/dPHI)(dU/dPHI)
    C
    C             + LAMBDA*U = F(THETA,PHI)
    C
    C     Where THETA is colatitude and PHI is longitude.
    C
    C     * * * * * * * *    Parameter Description     * * * * * * * * * *
    C
    C             * * * * * *   On Input    * * * * * *
    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     * * * * * * * * * * * * * * 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
    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       M must be greater than 5.
    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     * * * * * * * * * * * * * * IMPORTANT * * * * * * * * * * * * * *
    C
    C     If PF is equal to 2*PI then it must be computed using the
    C     statement PF = 2.*PIMACH(DUM). This insures that PF in the users
    C     program is equal to 2*PI in this program which permits tests of
    C     the input parameters that otherwise would not be possible.
    C
    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.
    C       N must be greater than 4.
    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,J) = U(I,N+J).
    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, HWSSSP 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
    C     IDIMF
    C       The row (or first) dimension of the array F as it appears in the
    C       program calling HWSSSP.  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. W may require up to 4*(N+1)+(16+INT(log2(N+1)))(M+1)
    C       locations. The actual number of locations used is computed by
    C       HWSSSP and is output in location W(1). INT( ) denotes the
    C       FORTRAN integer function.
    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.  HWSSSP 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.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 a possibly incorrect
    C       call to HWSSSP, the user should test IERROR after a call.
    C
    C     W
    C       Contains intermediate values that must not be destroyed if
    C       HWSSSP will be called again with INTL = 1. W(1) contains the
    C       required length of W .
    C
    C *Long Description:
    C
    C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
    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         January 1978
    C     Revision
    C
    C
    C     Subprograms    HWSSSP,HWSSS1,GENBUN,POISD2,POISN2,POISP2,COSGEN,ME
    C     Required       TRIX,TRI3,PIMACH
    C
    C     Special        NONE
    C     Conditions
    C
    C     Common         NONE
    C     Blocks
    C
    C     I/O            NONE
    C
    C     Precision      Single
    C
    C     Specialist     Paul Swarztrauber
    C
    C     Language       FORTRAN
    C
    C     History        Version 1 - September 1973
    C                    Version 2 - April     1976
    C                    Version 3 - January   1978
    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                    GENBUN to solve the system.
    C
    C     Space
    C     Required       CONTROL DATA 7600
    C
    C     Timing and        The execution time T on the NCAR Control Data
    C     Accuracy       7600 for subroutine HWSSSP is roughly proportional
    C                    to M*N*log2(N), but also depends on the input
    C                    parameters NBDCND and MBDCND.  Some typical values
    C                    are listed in the table below.
    C                       The solution process employed results in a loss
    C                    of no more than three significant digits for N and
    C                    M as large as 64.  More detailed information 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
    C                       M(=N)    MBDCND    NBDCND    T(MSECS)
    C                       -----    ------    ------    --------
    C
    C                        32        0         0          31
    C                        32        1         1          23
    C                        32        3         3          36
    C                        64        0         0         128
    C                        64        1         1          96
    C                        64        3         3         142
    C
    C     Portability    American National Standards Institute FORTRAN.
    C                    The machine dependent constant PI is defined in
    C                    function PIMACH.
    C
    C     Required       SIN,COS
    C     Resident
    C     Routines
    C
    C     References     P. N. Swarztrauber,'The Direct Solution Of The
    C                    Discrete Poisson Equation On The Surface Of a
    C                    Sphere, SIAM J. Numer. Anal.,15(1974), pp 212-215
    C
    C                    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               P. N. Swarztrauber, The direct solution of the discrete
    C                 Poisson equation on the surface of a sphere, SIAM
    C                 Journal on Numerical Analysis 15 (1974), pp. 212-215.
    C***ROUTINES CALLED  HWSSS1, PIMACH
    C***REVISION HISTORY  (YYMMDD)
    C   801001  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   891009  Removed unreferenced variable.  (WRB)
    C   891009  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  HWSSSP