*DECK HSTPLR
      SUBROUTINE HSTPLR (A, B, M, MBDCND, BDA, BDB, C, D, N, NBDCND,
     +   BDC, BDD, ELMBDA, F, IDIMF, PERTRB, IERROR, W)
C***BEGIN PROLOGUE  HSTPLR
C***PURPOSE  Solve the standard five-point finite difference
C            approximation on a staggered grid to the Helmholtz equation
C            in polar coordinates.
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B1A1A
C***TYPE      SINGLE PRECISION (HSTPLR-S)
C***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, POLAR
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C      HSTPLR solves the standard five-point finite difference
C      approximation on a staggered grid to the Helmholtz equation in
C      polar coordinates
C
C      (1/R)(d/DR)(R(dU/DR)) + (1/R**2)(d/dTHETA)(dU/dTHETA)
C
C                      + LAMBDA*U = F(R,THETA)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C             * * * * * *   On Input    * * * * * *
C
C    A,B
C      The range of R, i.e. A .LE. R .LE. B.  A must be less than B and
C      A must be non-negative.
C
C    M
C      The number of grid points in the interval (A,B).  The grid points
C      in the R-direction are given by R(I) = A + (I-0.5)DR for
C      I=1,2,...,M where DR =(B-A)/M.  M must be greater than 2.
C
C    MBDCND
C      Indicates the type of boundary conditions at R = A and R = B.
C
C      = 1  If the solution is specified at R = A and R = B.
C
C      = 2  If the solution is specified at R = A and the derivative
C           of the solution with respect to R is specified at R = B.
C           (see note 1 below)
C
C      = 3  If the derivative of the solution with respect to R is
C           specified at R = A (see note 2 below) and R = B.
C
C      = 4  If the derivative of the solution with respect to R is
C           specified at R = A (see note 2 below) and the solution is
C           specified at R = B.
C
C      = 5  If the solution is unspecified at R = A = 0 and the solution
C           is specified at R = B.
C
C      = 6  If the solution is unspecified at R = A = 0 and the
C           derivative of the solution with respect to R is specified at
C           R = B.
C
C      NOTE 1:  If A = 0, MBDCND = 2, and NBDCND = 0 or 3, the system of
C               equations to be solved is singular.  The unique solution
C               is determined by extrapolation to the specification of
C               U(0,THETA(1)).  But in this case the right side of the
C               system will be perturbed by the constant PERTRB.
C
C      NOTE 2:  If A = 0, do not use MBDCND = 3 or 4, but instead use
C               MBDCND = 1,2,5, or 6.
C
C    BDA
C      A one-dimensional array of length N that specifies the boundary
C      values (if any) of the solution at R = A.  When MBDCND = 1 or 2,
C
C               BDA(J) = U(A,THETA(J)) ,          J=1,2,...,N.
C
C      When MBDCND = 3 or 4,
C
C               BDA(J) = (d/dR)U(A,THETA(J)) ,    J=1,2,...,N.
C
C      When MBDCND = 5 or 6, BDA is a dummy variable.
C
C    BDB
C      A one-dimensional array of length N that specifies the boundary
C      values of the solution at R = B.  When MBDCND = 1,4, or 5,
C
C               BDB(J) = U(B,THETA(J)) ,          J=1,2,...,N.
C
C      When MBDCND = 2,3, or 6,
C
C               BDB(J) = (d/dR)U(B,THETA(J)) ,    J=1,2,...,N.
C
C    C,D
C      The range of THETA, i.e. C .LE. THETA .LE. D.  C must be less
C      than D.
C
C    N
C      The number of unknowns in the interval (C,D).  The unknowns in
C      the THETA-direction are given by THETA(J) = C + (J-0.5)DT,
C      J=1,2,...,N, where DT = (D-C)/N.  N must be greater than 2.
C
C    NBDCND
C      Indicates the type of boundary conditions at THETA = C
C      and THETA = D.
C
C      = 0  If the solution is periodic in THETA, i.e.
C           U(I,J) = U(I,N+J).
C
C      = 1  If the solution is specified at THETA = C and THETA = D
C           (see note below).
C
C      = 2  If the solution is specified at THETA = C and the derivative
C           of the solution with respect to THETA is specified at
C           THETA = D (see note below).
C
C      = 3  If the derivative of the solution with respect to THETA is
C           specified at THETA = C and THETA = D.
C
C      = 4  If the derivative of the solution with respect to THETA is
C           specified at THETA = C and the solution is specified at
C           THETA = d (see note below).
C
C      NOTE:  When NBDCND = 1, 2, or 4, do not use MBDCND = 5 or 6 (the
C      former indicates that the solution is specified at R =  0; the
C      latter indicates the solution is unspecified at R = 0).  Use
C      instead MBDCND = 1 or 2.
C
C    BDC
C      A one dimensional array of length M that specifies the boundary
C      values of the solution at THETA = C.   When NBDCND = 1 or 2,
C
C               BDC(I) = U(R(I),C) ,              I=1,2,...,M.
C
C      When NBDCND = 3 or 4,
C
C               BDC(I) = (d/dTHETA)U(R(I),C),     I=1,2,...,M.
C
C      When NBDCND = 0, BDC is a dummy variable.
C
C    BDD
C      A one-dimensional array of length M that specifies the boundary
C      values of the solution at THETA = D.  When NBDCND = 1 or 4,
C
C               BDD(I) = U(R(I),D) ,              I=1,2,...,M.
C
C      When NBDCND = 2 or 3,
C
C               BDD(I) = (d/dTHETA)U(R(I),D) ,    I=1,2,...,M.
C
C      When NBDCND = 0, BDD is a dummy variable.
C
C    ELMBDA
C      The constant LAMBDA in the Helmholtz equation.  If LAMBDA is
C      greater than 0, a solution may not exist.  However, HSTPLR will
C      attempt to find a solution.
C
C    F
C      A two-dimensional array that specifies the values of the right
C      side of the Helmholtz equation.  For I=1,2,...,M and J=1,2,...,N
C
C               F(I,J) = F(R(I),THETA(J)) .
C
C      F must be dimensioned at least M X N.
C
C    IDIMF
C      The row (or first) dimension of the array F as it appears in the
C      program calling HSTPLR.  This parameter is used to specify the
C      variable dimension of F.  IDIMF must be at least M.
C
C    W
C      A one-dimensional array that must be provided by the user for
C      work space.  W may require up to 13M + 4N + M*INT(log2(N))
C      locations.  The actual number of locations used is computed by
C      HSTPLR and is returned in the location W(1).
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 (R(I),THETA(J)) for
C      I=1,2,...,M, J=1,2,...,N.
C
C    PERTRB
C      If a combination of periodic, derivative, or unspecified
C      boundary conditions is specified for a Poisson equation
C      (LAMBDA = 0), a solution may not exist.  PERTRB is a con-
C      stant, calculated and subtracted from F, which ensures
C      that a solution exists.  HSTPLR then computes this
C      solution, which is a least squares solution to the
C      original approximation.  This solution plus any constant is also
C      a solution; hence, the solution is not unique.  The value of
C      PERTRB should be small compared to the right side F.
C      Otherwise, a solution is obtained to an essentially different
C      problem.  This comparison should always be made to insure that
C      a meaningful solution has been obtained.
C
C    IERROR
C      An error flag that indicates invalid input parameters.
C      Except for numbers 0 and 11, a solution is not attempted.
C
C      =  0  No error
C
C      =  1  A .LT. 0
C
C      =  2  A .GE. B
C
C      =  3  MBDCND .LT. 1 or MBDCND .GT. 6
C
C      =  4  C .GE. D
C
C      =  5  N .LE. 2
C
C      =  6  NBDCND .LT. 0 or NBDCND .GT. 4
C
C      =  7  A = 0 and MBDCND = 3 or 4
C
C      =  8  A .GT. 0 and MBDCND .GE. 5
C
C      =  9  MBDCND .GE. 5 and NBDCND .NE. 0 or 3
C
C      = 10  IDIMF .LT. M
C
C      = 11  LAMBDA .GT. 0
C
C      = 12  M .LE. 2
C
C      Since this is the only means of indicating a possibly
C      incorrect call to HSTPLR, the user should test IERROR after
C      the call.
C
C    W
C      W(1) contains the required length of W.
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
C     Arguments      W(see ARGUMENT LIST)
C
C     Latest         June 1, 1977
C     Revision
C
C     Subprograms    HSTPLR,POISTG,POSTG2,GENBUN,POISD2,POISN2,POISP2,
C     Required       COSGEN,MERGE,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     Roland Sweet
C
C     Language       FORTRAN
C
C     History        Written by Roland Sweet at NCAR in February, 1977
C
C     Algorithm      This subroutine defines the finite-difference
C                    equations, incorporates boundary data, adjusts the
C                    right side when the system is singular and calls
C                    either POISTG or GENBUN which solves the linear
C                    system of equations.
C
C     Space          8265(decimal) = 20111(octal) LOCATIONS ON THE
C     Required       NCAR Control Data 7600
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine HSTPLR is roughly proportional
C                    to M*N*log2(N).  Some typical values are listed in
C                    the table below.
C                       The solution process employed results in a loss
C                    of no more than four significant digits for N and M
C                    as large as 64.  More detailed information about
C                    accuracy can be found in the documentation for
C                    subroutine POISTG which is the routine that
C                    actually solves the finite difference equations.
C
C
C                       M(=N)    MBDCND    NBDCND    T(MSECS)
C                       -----    ------    ------    --------
C
C                        32       1-6       1-4         56
C                        64       1-6       1-4        230
C
C     Portability    American National Standards Institute Fortran.
C                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       COS
C     Resident
C     Routines
C
C     Reference      Schumann, U. 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
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES  U. Schumann and R. Sweet, A direct method for the
C                 solution of Poisson's equation with Neumann boundary
C                 conditions on a staggered grid of arbitrary size,
C                 Journal of Computational Physics 20, (1976),
C                 pp. 171-182.
C***ROUTINES CALLED  GENBUN, POISTG
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  HSTPLR
C
C
      DIMENSION       F(IDIMF,*)
      DIMENSION       BDA(*)     ,BDB(*)     ,BDC(*)     ,BDD(*)     ,
     1                W(*)
C***FIRST EXECUTABLE STATEMENT  HSTPLR
      IERROR = 0
      IF (A .LT. 0.) IERROR = 1
      IF (A .GE. B) IERROR = 2
      IF (MBDCND.LE.0 .OR. MBDCND.GE.7) IERROR = 3
      IF (C .GE. D) IERROR = 4
      IF (N .LE. 2) IERROR = 5
      IF (NBDCND.LT.0 .OR. NBDCND.GE.5) IERROR = 6
      IF (A.EQ.0. .AND. (MBDCND.EQ.3 .OR. MBDCND.EQ.4)) IERROR = 7
      IF (A.GT.0. .AND. MBDCND.GE.5) IERROR = 8
      IF (MBDCND.GE.5 .AND. NBDCND.NE.0 .AND. NBDCND.NE.3) IERROR = 9
      IF (IDIMF .LT. M) IERROR = 10
      IF (M .LE. 2) IERROR = 12
      IF (IERROR .NE. 0) RETURN
      DELTAR = (B-A)/M
      DLRSQ = DELTAR**2
      DELTHT = (D-C)/N
      DLTHSQ = DELTHT**2
      NP = NBDCND+1
      ISW = 1
      MB = MBDCND
      IF (A.EQ.0. .AND. MBDCND.EQ.2) MB = 6
C
C     DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
C
      IWB = M
      IWC = IWB+M
      IWR = IWC+M
      DO 101 I=1,M
         J = IWR+I
         W(J) = A+(I-0.5)*DELTAR
         W(I) = (A+(I-1)*DELTAR)/DLRSQ
         K = IWC+I
         W(K) = (A+I*DELTAR)/DLRSQ
         K = IWB+I
         W(K) = (ELMBDA-2./DLRSQ)*W(J)
  101 CONTINUE
      DO 103 I=1,M
         J = IWR+I
         A1 = W(J)
         DO 102 J=1,N
            F(I,J) = A1*F(I,J)
  102    CONTINUE
  103 CONTINUE
C
C     ENTER BOUNDARY DATA FOR R-BOUNDARIES.
C
      GO TO (104,104,106,106,108,108),MB
  104 A1 = 2.*W(1)
      W(IWB+1) = W(IWB+1)-W(1)
      DO 105 J=1,N
         F(1,J) = F(1,J)-A1*BDA(J)
  105 CONTINUE
      GO TO 108
  106 A1 = DELTAR*W(1)
      W(IWB+1) = W(IWB+1)+W(1)
      DO 107 J=1,N
         F(1,J) = F(1,J)+A1*BDA(J)
  107 CONTINUE
  108 GO TO (109,111,111,109,109,111),MB
  109 A1 = 2.*W(IWR)
      W(IWC) = W(IWC)-W(IWR)
      DO 110 J=1,N
         F(M,J) = F(M,J)-A1*BDB(J)
  110 CONTINUE
      GO TO 113
  111 A1 = DELTAR*W(IWR)
      W(IWC) = W(IWC)+W(IWR)
      DO 112 J=1,N
         F(M,J) = F(M,J)-A1*BDB(J)
  112 CONTINUE
C
C     ENTER BOUNDARY DATA FOR THETA-BOUNDARIES.
C
  113 A1 = 2./DLTHSQ
      GO TO (123,114,114,116,116),NP
  114 DO 115 I=1,M
         J = IWR+I
         F(I,1) = F(I,1)-A1*BDC(I)/W(J)
  115 CONTINUE
      GO TO 118
  116 A1 = 1./DELTHT
      DO 117 I=1,M
         J = IWR+I
         F(I,1) = F(I,1)+A1*BDC(I)/W(J)
  117 CONTINUE
  118 A1 = 2./DLTHSQ
      GO TO (123,119,121,121,119),NP
  119 DO 120 I=1,M
         J = IWR+I
         F(I,N) = F(I,N)-A1*BDD(I)/W(J)
  120 CONTINUE
      GO TO 123
  121 A1 = 1./DELTHT
      DO 122 I=1,M
         J = IWR+I
         F(I,N) = F(I,N)-A1*BDD(I)/W(J)
  122 CONTINUE
  123 CONTINUE
C
C     ADJUST RIGHT SIDE OF SINGULAR PROBLEMS TO INSURE EXISTENCE OF A
C     SOLUTION.
C
      PERTRB = 0.
      IF (ELMBDA) 133,125,124
  124 IERROR = 11
      GO TO 133
  125 GO TO (133,133,126,133,133,126),MB
  126 GO TO (127,133,133,127,133),NP
  127 CONTINUE
      ISW = 2
      DO 129 J=1,N
         DO 128 I=1,M
            PERTRB = PERTRB+F(I,J)
  128    CONTINUE
  129 CONTINUE
      PERTRB = PERTRB/(M*N*0.5*(A+B))
      DO 131 I=1,M
         J = IWR+I
         A1 = PERTRB*W(J)
         DO 130 J=1,N
            F(I,J) = F(I,J)-A1
  130    CONTINUE
  131 CONTINUE
      A2 = 0.
      DO 132 J=1,N
         A2 = A2+F(1,J)
  132 CONTINUE
      A2 = A2/W(IWR+1)
  133 CONTINUE
C
C     MULTIPLY I-TH EQUATION THROUGH BY  R(I)*DELTHT**2
C
      DO 135 I=1,M
         J = IWR+I
         A1 = DLTHSQ*W(J)
         W(I) = A1*W(I)
         J = IWC+I
         W(J) = A1*W(J)
         J = IWB+I
         W(J) = A1*W(J)
         DO 134 J=1,N
            F(I,J) = A1*F(I,J)
  134    CONTINUE
  135 CONTINUE
      LP = NBDCND
      W(1) = 0.
      W(IWR) = 0.
C
C     CALL POISTG OR GENBUN TO SOLVE THE SYSTEM OF EQUATIONS.
C
      IF (LP .EQ. 0) GO TO 136
      CALL POISTG (LP,N,1,M,W,W(IWB+1),W(IWC+1),IDIMF,F,IERR1,W(IWR+1))
      GO TO 137
  136 CALL GENBUN (LP,N,1,M,W,W(IWB+1),W(IWC+1),IDIMF,F,IERR1,W(IWR+1))
  137 CONTINUE
      W(1) = W(IWR+1)+3*M
      IF (A.NE.0. .OR. MBDCND.NE.2 .OR. ISW.NE.2) GO TO 141
      A1 = 0.
      DO 138 J=1,N
         A1 = A1+F(1,J)
  138 CONTINUE
      A1 = (A1-DLRSQ*A2/16.)/N
      IF (NBDCND .EQ. 3) A1 = A1+(BDD(1)-BDC(1))/(D-C)
      A1 = BDA(1)-A1
      DO 140 I=1,M
         DO 139 J=1,N
            F(I,J) = F(I,J)+A1
  139    CONTINUE
  140 CONTINUE
  141 CONTINUE
      RETURN
      END
*DECK POISTG
      SUBROUTINE POISTG (NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y,
     +   IERROR, W)
C***BEGIN PROLOGUE  POISTG
C***PURPOSE  Solve a block tridiagonal system of linear equations
C            that results from a staggered grid finite difference
C            approximation to 2-D elliptic PDE's.
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B4B
C***TYPE      SINGLE PRECISION (POISTG-S)
C***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE, TRIDIAGONAL
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine POISTG solves the linear system of equations
C
C       A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
C       + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
C
C       for I=1,2,...,M and J=1,2,...,N.
C
C     The indices I+1 and I-1 are evaluated modulo M, i.e.
C     X(0,J) = X(M,J) and X(M+1,J) = X(1,J), and X(I,0) may be equal to
C     X(I,1) or -X(I,1) and X(I,N+1) may be equal to X(I,N) or -X(I,N)
C     depending on an input parameter.
C
C
C     * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C             * * * * * *   On Input    * * * * * *
C
C   NPEROD
C     Indicates the values which X(I,0) and X(I,N+1) are assumed
C     to have.
C     = 1 If X(I,0) = -X(I,1) and X(I,N+1) = -X(I,N)
C     = 2 If X(I,0) = -X(I,1) and X(I,N+1) =  X(I,N)
C     = 3 If X(I,0) =  X(I,1) and X(I,N+1) =  X(I,N)
C     = 4 If X(I,0) =  X(I,1) and X(I,N+1) = -X(I,N)
C
C   N
C     The number of unknowns in the J-direction.  N must
C     be greater than 2.
C
C   MPEROD
C     = 0 If A(1) and C(M) are not zero
C     = 1 If A(1) = C(M) = 0
C
C   M
C     The number of unknowns in the I-direction.  M must
C     be greater than 2.
C
C   A,B,C
C     One-dimensional arrays of length M that specify the coefficients
C     in the linear equations given above.  If MPEROD = 0 the array
C     elements must not depend on the index I, but must be constant.
C     Specifically, the subroutine checks the following condition
C
C           A(I) = C(1)
C           B(I) = B(1)
C           C(I) = C(1)
C
C     for I = 1, 2, ..., M.
C
C   IDIMY
C     The row (or first) dimension of the two-dimensional array Y as
C     it appears in the program calling POISTG.  This parameter is
C     used to specify the variable dimension of Y.  IDIMY must be at
C     least M.
C
C   Y
C     A two-dimensional array that specifies the values of the
C     right side of the linear system of equations given above.
C     Y must be dimensioned at least M X N.
C
C   W
C     A one-dimensional work array that must be provided by the user
C     for work space.  W may require up to 9M + 4N + M(INT(log2(N)))
C     locations.  The actual number of locations used is computed by
C     POISTG and returned in location W(1).
C
C
C             * * * * * *   On Output     * * * * * *
C
C   Y
C     Contains the solution X.
C
C   IERROR
C     An error flag that indicates invalid input parameters.  Except
C     for number zero, a solution is not attempted.
C     = 0  No error
C     = 1  If M .LE. 2
C     = 2  If N .LE. 2
C     = 3  IDIMY .LT. M
C     = 4  If NPEROD .LT. 1 or NPEROD .GT. 4
C     = 5  If MPEROD .LT. 0 or MPEROD .GT. 1
C     = 6  If MPEROD = 0 and
C          A(I) .NE. C(1) or B(I) .NE. B(1) or C(I) .NE. C(1)
C          for some I = 1, 2, ..., M.
C       = 7 If MPEROD .EQ. 1 .AND. (A(1).NE.0 .OR. C(M).NE.0)
C
C   W
C     W(1) contains the required length of W.
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   A(M),B(M),C(M),Y(IDIMY,N),
C     Arguments      W(see argument list)
C
C     Latest         June 1, 1977
C     Revision
C
C     Subprograms    POISTG,POSTG2,COSGEN,MERGE,TRIX,TRI3,PIMACH
C     Required
C
C     Special        NONE
C     Conditions
C
C     Common         NONE
C     Blocks
C
C     I/O            NONE
C
C     Precision      Single
C
C     Specialist     Roland Sweet
C
C     Language       FORTRAN
C
C     History        Written by Roland Sweet in 1973
C                    Revised by Roland Sweet in 1977
C
C
C     Space          3297(decimal) = 6341(octal) locations on the
C     Required       NCAR Control Data 7600
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine POISTG is roughly proportional
C                    to M*N*log2(N).  Some typical values are listed
C                    in the table below.  More comprehensive timing
C                    charts may be found in the reference.
C                       To measure the accuracy of the algorithm a
C                    uniform random number generator was used to create
C                    a solution array X for the system given in the
C                    'PURPOSE ' with
C
C                       A(I) = C(I) = -0.5*B(I) = 1,       I=1,2,...,M
C
C                    and, when MPEROD = 1
C
C                       A(1) = C(M) = 0
C                       B(1) = B(M) =-1.
C
C                    The solution X was substituted into the given sys-
C                    tem and, using double precision, a right side Y was
C                    computed.  Using this array Y subroutine POISTG was
C                    called to produce an approximate solution Z.  Then
C                    the relative error, defined as
C
C                       E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J)))
C
C                    where the two maxima are taken over all I=1,2,...,M
C                    and J=1,2,...,N, was computed.  The value of E is
C                    given in the table below for some typical values of
C                    M and N.
C
C
C                       M (=N)    MPEROD    NPEROD    T(MSECS)    E
C                       ------    ------    ------    --------  ------
C
C                         31        0-1       1-4        45     9.E-13
C                         31        1         1          21     4.E-13
C                         31        1         3          41     3.E-13
C                         32        0-1       1-4        51     3.E-12
C                         32        1         1          32     3.E-13
C                         32        1         3          48     1.E-13
C                         33        0-1       1-4        42     1.E-12
C                         33        1         1          30     4.E-13
C                         33        1         3          34     1.E-13
C                         63        0-1       1-4       186     3.E-12
C                         63        1         1          91     1.E-12
C                         63        1         3         173     2.E-13
C                         64        0-1       1-4       209     4.E-12
C                         64        1         1         128     1.E-12
C                         64        1         3         199     6.E-13
C                         65        0-1       1-4       143     2.E-13
C                         65        1         1         160     1.E-11
C                         65        1         3         138     4.E-13
C
C     Portability    American National Standards Institute FORTRAN.
C                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       COS
C     Resident
C     Routines
C
C     Reference      Schumann, U. 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
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES  U. Schumann and R. Sweet, A direct method for the
C                 solution of Poisson's equation with Neumann boundary
C                 conditions on a staggered grid of arbitrary size,
C                 Journal of Computational Physics 20, (1976),
C                 pp. 171-182.
C***ROUTINES CALLED  POSTG2
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   861211  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  POISTG
C
C
      DIMENSION       Y(IDIMY,*)
      DIMENSION       W(*)       ,B(*)       ,A(*)       ,C(*)
C***FIRST EXECUTABLE STATEMENT  POISTG
      IERROR = 0
      IF (M .LE. 2) IERROR = 1
      IF (N .LE. 2) IERROR = 2
      IF (IDIMY .LT. M) IERROR = 3
      IF (NPEROD.LT.1 .OR. NPEROD.GT.4) IERROR = 4
      IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR = 5
      IF (MPEROD .EQ. 1) GO TO 103
      DO 101 I=1,M
         IF (A(I) .NE. C(1)) GO TO 102
         IF (C(I) .NE. C(1)) GO TO 102
         IF (B(I) .NE. B(1)) GO TO 102
  101 CONTINUE
      GO TO 104
  102 IERROR = 6
      RETURN
  103 IF (A(1).NE.0. .OR. C(M).NE.0.) IERROR = 7
  104 IF (IERROR .NE. 0) RETURN
      IWBA = M+1
      IWBB = IWBA+M
      IWBC = IWBB+M
      IWB2 = IWBC+M
      IWB3 = IWB2+M
      IWW1 = IWB3+M
      IWW2 = IWW1+M
      IWW3 = IWW2+M
      IWD = IWW3+M
      IWTCOS = IWD+M
      IWP = IWTCOS+4*N
      DO 106 I=1,M
         K = IWBA+I-1
         W(K) = -A(I)
         K = IWBC+I-1
         W(K) = -C(I)
         K = IWBB+I-1
         W(K) = 2.-B(I)
         DO 105 J=1,N
            Y(I,J) = -Y(I,J)
  105    CONTINUE
  106 CONTINUE
      NP = NPEROD
      MP = MPEROD+1
      GO TO (110,107),MP
  107 CONTINUE
      GO TO (108,108,108,119),NPEROD
  108 CONTINUE
      CALL POSTG2 (NP,N,M,W(IWBA),W(IWBB),W(IWBC),IDIMY,Y,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
      IPSTOR = W(IWW1)
      IREV = 2
      IF (NPEROD .EQ. 4) GO TO 120
  109 CONTINUE
      GO TO (123,129),MP
  110 CONTINUE
C
C     REORDER UNKNOWNS WHEN MP =0
C
      MH = (M+1)/2
      MHM1 = MH-1
      MODD = 1
      IF (MH*2 .EQ. M) MODD = 2
      DO 115 J=1,N
         DO 111 I=1,MHM1
            MHPI = MH+I
            MHMI = MH-I
            W(I) = Y(MHMI,J)-Y(MHPI,J)
            W(MHPI) = Y(MHMI,J)+Y(MHPI,J)
  111    CONTINUE
         W(MH) = 2.*Y(MH,J)
         GO TO (113,112),MODD
  112    W(M) = 2.*Y(M,J)
  113    CONTINUE
         DO 114 I=1,M
            Y(I,J) = W(I)
  114    CONTINUE
  115 CONTINUE
      K = IWBC+MHM1-1
      I = IWBA+MHM1
      W(K) = 0.
      W(I) = 0.
      W(K+1) = 2.*W(K+1)
      GO TO (116,117),MODD
  116 CONTINUE
      K = IWBB+MHM1-1
      W(K) = W(K)-W(I-1)
      W(IWBC-1) = W(IWBC-1)+W(IWBB-1)
      GO TO 118
  117 W(IWBB-1) = W(K+1)
  118 CONTINUE
      GO TO 107
  119 CONTINUE
C
C     REVERSE COLUMNS WHEN NPEROD = 4.
C
      IREV = 1
      NBY2 = N/2
      NP = 2
  120 DO 122 J=1,NBY2
         MSKIP = N+1-J
         DO 121 I=1,M
            A1 = Y(I,J)
            Y(I,J) = Y(I,MSKIP)
            Y(I,MSKIP) = A1
  121    CONTINUE
  122 CONTINUE
      GO TO (108,109),IREV
  123 CONTINUE
      DO 128 J=1,N
         DO 124 I=1,MHM1
            MHMI = MH-I
            MHPI = MH+I
            W(MHMI) = .5*(Y(MHPI,J)+Y(I,J))
            W(MHPI) = .5*(Y(MHPI,J)-Y(I,J))
  124    CONTINUE
         W(MH) = .5*Y(MH,J)
         GO TO (126,125),MODD
  125    W(M) = .5*Y(M,J)
  126    CONTINUE
         DO 127 I=1,M
            Y(I,J) = W(I)
  127    CONTINUE
  128 CONTINUE
  129 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR W ARRAY.
C
      W(1) = IPSTOR+IWP-1
      RETURN
      END
*DECK POSTG2
      SUBROUTINE POSTG2 (NPEROD, N, M, A, BB, C, IDIMQ, Q, B, B2, B3, W,
     +   W2, W3, D, TCOS, P)
C***BEGIN PROLOGUE  POSTG2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to POISTG
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (POSTG2-S)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve Poisson's equation on a staggered grid.
C
C***SEE ALSO  POISTG
C***ROUTINES CALLED  COSGEN, S1MERG, TRI3, TRIX
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   920130  Modified to use merge routine S1MERG rather than deleted
C           routine MERGE.  (WRB)
C***END PROLOGUE  POSTG2
C
      DIMENSION       A(*)       ,BB(*)      ,C(*)       ,Q(IDIMQ,*) ,
     1                B(*)       ,B2(*)      ,B3(*)      ,W(*)       ,
     2                W2(*)      ,W3(*)      ,D(*)       ,TCOS(*)    ,
     3                K(4)       ,P(*)
      EQUIVALENCE     (K(1),K1)  ,(K(2),K2)  ,(K(3),K3)  ,(K(4),K4)
C***FIRST EXECUTABLE STATEMENT  POSTG2
      NP = NPEROD
      FNUM = 0.5*(NP/3)
      FNUM2 = 0.5*(NP/2)
      MR = M
      IP = -MR
      IPSTOR = 0
      I2R = 1
      JR = 2
      NR = N
      NLAST = N
      KR = 1
      LR = 0
      IF (NR .LE. 3) GO TO 142
  101 CONTINUE
      JR = 2*I2R
      NROD = 1
      IF ((NR/2)*2 .EQ. NR) NROD = 0
      JSTART = 1
      JSTOP = NLAST-JR
      IF (NROD .EQ. 0) JSTOP = JSTOP-I2R
      I2RBY2 = I2R/2
      IF (JSTOP .GE. JSTART) GO TO 102
      J = JR
      GO TO 115
  102 CONTINUE
C
C     REGULAR REDUCTION.
C
      IJUMP = 1
      DO 114 J=JSTART,JSTOP,JR
         JP1 = J+I2RBY2
         JP2 = J+I2R
         JP3 = JP2+I2RBY2
         JM1 = J-I2RBY2
         JM2 = J-I2R
         JM3 = JM2-I2RBY2
         IF (J .NE. 1) GO TO 106
         CALL COSGEN (I2R,1,FNUM,0.5,TCOS)
         IF (I2R .NE. 1) GO TO 104
         DO 103 I=1,MR
            B(I) = Q(I,1)
            Q(I,1) = Q(I,2)
  103    CONTINUE
         GO TO 112
  104    DO 105 I=1,MR
            B(I) = Q(I,1)+0.5*(Q(I,JP2)-Q(I,JP1)-Q(I,JP3))
            Q(I,1) = Q(I,JP2)+Q(I,1)-Q(I,JP1)
  105    CONTINUE
         GO TO 112
  106    CONTINUE
         GO TO (107,108),IJUMP
  107    CONTINUE
         IJUMP = 2
         CALL COSGEN (I2R,1,0.5,0.0,TCOS)
  108    CONTINUE
         IF (I2R .NE. 1) GO TO 110
         DO 109 I=1,MR
            B(I) = 2.*Q(I,J)
            Q(I,J) = Q(I,JM2)+Q(I,JP2)
  109    CONTINUE
         GO TO 112
  110    DO 111 I=1,MR
            FI = Q(I,J)
            Q(I,J) = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2)
            B(I) = FI+Q(I,J)-Q(I,JM3)-Q(I,JP3)
  111    CONTINUE
  112    CONTINUE
         CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
         DO 113 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  113    CONTINUE
C
C     END OF REDUCTION FOR REGULAR UNKNOWNS.
C
  114 CONTINUE
C
C     BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN.
C
      J = JSTOP+JR
  115 NLAST = J
      JM1 = J-I2RBY2
      JM2 = J-I2R
      JM3 = JM2-I2RBY2
      IF (NROD .EQ. 0) GO TO 125
C
C     ODD NUMBER OF UNKNOWNS
C
      IF (I2R .NE. 1) GO TO 117
      DO 116 I=1,MR
         B(I) = Q(I,J)
         Q(I,J) = Q(I,JM2)
  116 CONTINUE
      GO TO 123
  117 DO 118 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  118 CONTINUE
      IF (NRODPR .NE. 0) GO TO 120
      DO 119 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)
  119 CONTINUE
      IP = IP-MR
      GO TO 122
  120 CONTINUE
      DO 121 I=1,MR
         Q(I,J) = Q(I,J)-Q(I,JM1)+Q(I,JM2)
  121 CONTINUE
  122 IF (LR .EQ. 0) GO TO 123
      CALL COSGEN (LR,1,FNUM2,0.5,TCOS(KR+1))
  123 CONTINUE
      CALL COSGEN (KR,1,FNUM2,0.5,TCOS)
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 124 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
  124 CONTINUE
      KR = KR+I2R
      GO TO 141
  125 CONTINUE
C
C     EVEN NUMBER OF UNKNOWNS
C
      JP1 = J+I2RBY2
      JP2 = J+I2R
      IF (I2R .NE. 1) GO TO 129
      DO 126 I=1,MR
         B(I) = Q(I,J)
  126 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      IP = 0
      IPSTOR = MR
      DO 127 I=1,MR
         P(I) = B(I)
         B(I) = B(I)+Q(I,N)
  127 CONTINUE
      TCOS(1) = -1.+2*(NP/2)
      TCOS(2) = 0.
      CALL TRIX (1,1,MR,A,BB,C,B,TCOS,D,W)
      DO 128 I=1,MR
         Q(I,J) = Q(I,JM2)+P(I)+B(I)
  128 CONTINUE
      GO TO 140
  129 CONTINUE
      DO 130 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  130 CONTINUE
      IF (NRODPR .NE. 0) GO TO 132
      DO 131 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  131 CONTINUE
      GO TO 134
  132 CONTINUE
      DO 133 I=1,MR
         B(I) = B(I)+Q(I,JP2)-Q(I,JP1)
  133 CONTINUE
  134 CONTINUE
      CALL COSGEN (I2R,1,0.5,0.0,TCOS)
      CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
      IP = IP+MR
      IPSTOR = MAX(IPSTOR,IP+MR)
      DO 135 I=1,MR
         II = IP+I
         P(II) = B(I)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         B(I) = P(II)+Q(I,JP2)
  135 CONTINUE
      IF (LR .EQ. 0) GO TO 136
      CALL COSGEN (LR,1,FNUM2,0.5,TCOS(I2R+1))
      CALL S1MERG (TCOS,0,I2R,I2R,LR,KR)
      GO TO 138
  136 DO 137 I=1,I2R
         II = KR+I
         TCOS(II) = TCOS(I)
  137 CONTINUE
  138 CALL COSGEN (KR,1,FNUM2,0.5,TCOS)
      CALL TRIX (KR,KR,MR,A,BB,C,B,TCOS,D,W)
      DO 139 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)+B(I)
  139 CONTINUE
  140 CONTINUE
      LR = KR
      KR = KR+JR
  141 CONTINUE
      NR = (NLAST-1)/JR+1
      IF (NR .LE. 3) GO TO 142
      I2R = JR
      NRODPR = NROD
      GO TO 101
  142 CONTINUE
C
C      BEGIN SOLUTION
C
      J = 1+JR
      JM1 = J-I2R
      JP1 = J+I2R
      JM2 = NLAST-I2R
      IF (NR .EQ. 2) GO TO 180
      IF (LR .NE. 0) GO TO 167
      IF (N .NE. 3) GO TO 156
C
C     CASE N = 3.
C
      GO TO (143,148,143),NP
  143 DO 144 I=1,MR
         B(I) = Q(I,2)
         B2(I) = Q(I,1)+Q(I,3)
         B3(I) = 0.
  144 CONTINUE
      GO TO (146,146,145),NP
  145 TCOS(1) = -1.
      TCOS(2) = 1.
      K1 = 1
      GO TO 147
  146 TCOS(1) = -2.
      TCOS(2) = 1.
      TCOS(3) = -1.
      K1 = 2
  147 K2 = 1
      K3 = 0
      K4 = 0
      GO TO 150
  148 DO 149 I=1,MR
         B(I) = Q(I,2)
         B2(I) = Q(I,3)
         B3(I) = Q(I,1)
  149 CONTINUE
      CALL COSGEN (3,1,0.5,0.0,TCOS)
      TCOS(4) = -1.
      TCOS(5) = 1.
      TCOS(6) = -1.
      TCOS(7) = 1.
      K1 = 3
      K2 = 2
      K3 = 1
      K4 = 1
  150 CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 151 I=1,MR
         B(I) = B(I)+B2(I)+B3(I)
  151 CONTINUE
      GO TO (153,153,152),NP
  152 TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
  153 DO 154 I=1,MR
         Q(I,2) = B(I)
         B(I) = Q(I,1)+B(I)
  154 CONTINUE
      TCOS(1) = -1.+4.*FNUM
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 155 I=1,MR
         Q(I,1) = B(I)
  155 CONTINUE
      JR = 1
      I2R = 0
      GO TO 188
C
C     CASE N = 2**P+1
C
  156 CONTINUE
      DO 157 I=1,MR
         B(I) = Q(I,J)+Q(I,1)-Q(I,JM1)+Q(I,NLAST)-Q(I,JM2)
  157 CONTINUE
      GO TO (158,160,158),NP
  158 DO 159 I=1,MR
         B2(I) = Q(I,1)+Q(I,NLAST)+Q(I,J)-Q(I,JM1)-Q(I,JP1)
         B3(I) = 0.
  159 CONTINUE
      K1 = NLAST-1
      K2 = NLAST+JR-1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(NLAST))
      TCOS(K2) = 2*NP-4
      CALL COSGEN (JR,1,0.5-FNUM,0.5,TCOS(K2+1))
      K3 = (3-NP)/2
      CALL S1MERG (TCOS,K1,JR-K3,K2-K3,JR+K3,0)
      K1 = K1-1+K3
      CALL COSGEN (JR,1,FNUM,0.5,TCOS(K1+1))
      K2 = JR
      K3 = 0
      K4 = 0
      GO TO 162
  160 DO 161 I=1,MR
         FI = (Q(I,J)-Q(I,JM1)-Q(I,JP1))/2.
         B2(I) = Q(I,1)+FI
         B3(I) = Q(I,NLAST)+FI
  161 CONTINUE
      K1 = NLAST+JR-1
      K2 = K1+JR-1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K1+1))
      CALL COSGEN (NLAST,1,0.5,0.0,TCOS(K2+1))
      CALL S1MERG (TCOS,K1,JR-1,K2,NLAST,0)
      K3 = K1+NLAST-1
      K4 = K3+JR
      CALL COSGEN (JR,1,0.5,0.5,TCOS(K3+1))
      CALL COSGEN (JR,1,0.0,0.5,TCOS(K4+1))
      CALL S1MERG (TCOS,K3,JR,K4,JR,K1)
      K2 = NLAST-1
      K3 = JR
      K4 = JR
  162 CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 163 I=1,MR
         B(I) = B(I)+B2(I)+B3(I)
  163 CONTINUE
      IF (NP .NE. 3) GO TO 164
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
  164 DO 165 I=1,MR
         Q(I,J) = B(I)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         B(I) = Q(I,J)+Q(I,1)
  165 CONTINUE
      CALL COSGEN (JR,1,FNUM,0.5,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      DO 166 I=1,MR
         Q(I,1) = Q(I,1)-Q(I,JM1)+B(I)
  166 CONTINUE
      GO TO 188
C
C     CASE OF GENERAL N WITH NR = 3 .
C
  167 CONTINUE
      DO 168 I=1,MR
         B(I) = Q(I,1)-Q(I,JM1)+Q(I,J)
  168 CONTINUE
      IF (NROD .NE. 0) GO TO 170
      DO 169 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  169 CONTINUE
      GO TO 172
  170 DO 171 I=1,MR
         B(I) = B(I)+Q(I,NLAST)-Q(I,JM2)
  171 CONTINUE
  172 CONTINUE
      DO 173 I=1,MR
         T = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         Q(I,J) = T
         B2(I) = Q(I,NLAST)+T
         B3(I) = Q(I,1)+T
  173 CONTINUE
      K1 = KR+2*JR
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K1+1))
      K2 = K1+JR
      TCOS(K2) = 2*NP-4
      K4 = (NP-1)*(3-NP)
      K3 = K2+1-K4
      CALL COSGEN (KR+JR+K4,1,K4/2.,1.-K4,TCOS(K3))
      K4 = 1-NP/3
      CALL S1MERG (TCOS,K1,JR-K4,K2-K4,KR+JR+K4,0)
      IF (NP .EQ. 3) K1 = K1-1
      K2 = KR+JR
      K4 = K1+K2
      CALL COSGEN (KR,1,FNUM2,0.5,TCOS(K4+1))
      K3 = K4+KR
      CALL COSGEN (JR,1,FNUM,0.5,TCOS(K3+1))
      CALL S1MERG (TCOS,K4,KR,K3,JR,K1)
      K4 = K3+JR
      CALL COSGEN (LR,1,FNUM2,0.5,TCOS(K4+1))
      CALL S1MERG (TCOS,K3,JR,K4,LR,K1+K2)
      CALL COSGEN (KR,1,FNUM2,0.5,TCOS(K3+1))
      K3 = KR
      K4 = KR
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 174 I=1,MR
         B(I) = B(I)+B2(I)+B3(I)
  174 CONTINUE
      IF (NP .NE. 3) GO TO 175
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
  175 DO 176 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
         B(I) = Q(I,1)+Q(I,J)
  176 CONTINUE
      CALL COSGEN (JR,1,FNUM,0.5,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      IF (JR .NE. 1) GO TO 178
      DO 177 I=1,MR
         Q(I,1) = B(I)
  177 CONTINUE
      GO TO 188
  178 CONTINUE
      DO 179 I=1,MR
         Q(I,1) = Q(I,1)-Q(I,JM1)+B(I)
  179 CONTINUE
      GO TO 188
  180 CONTINUE
C
C     CASE OF GENERAL N AND NR = 2 .
C
      DO 181 I=1,MR
         II = IP+I
         B3(I) = 0.
         B(I) = Q(I,1)+P(II)
         Q(I,1) = Q(I,1)-Q(I,JM1)
         B2(I) = Q(I,1)+Q(I,NLAST)
  181 CONTINUE
      K1 = KR+JR
      K2 = K1+JR
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K1+1))
      GO TO (182,183,182),NP
  182 TCOS(K2) = 2*NP-4
      CALL COSGEN (KR,1,0.0,1.0,TCOS(K2+1))
      GO TO 184
  183 CALL COSGEN (KR+1,1,0.5,0.0,TCOS(K2))
  184 K4 = 1-NP/3
      CALL S1MERG (TCOS,K1,JR-K4,K2-K4,KR+K4,0)
      IF (NP .EQ. 3) K1 = K1-1
      K2 = KR
      CALL COSGEN (KR,1,FNUM2,0.5,TCOS(K1+1))
      K4 = K1+KR
      CALL COSGEN (LR,1,FNUM2,0.5,TCOS(K4+1))
      K3 = LR
      K4 = 0
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 185 I=1,MR
         B(I) = B(I)+B2(I)
  185 CONTINUE
      IF (NP .NE. 3) GO TO 186
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
  186 DO 187 I=1,MR
         Q(I,1) = Q(I,1)+B(I)
  187 CONTINUE
  188 CONTINUE
C
C     START BACK SUBSTITUTION.
C
      J = NLAST-JR
      DO 189 I=1,MR
         B(I) = Q(I,NLAST)+Q(I,J)
  189 CONTINUE
      JM2 = NLAST-I2R
      IF (JR .NE. 1) GO TO 191
      DO 190 I=1,MR
         Q(I,NLAST) = 0.
  190 CONTINUE
      GO TO 195
  191 CONTINUE
      IF (NROD .NE. 0) GO TO 193
      DO 192 I=1,MR
         II = IP+I
         Q(I,NLAST) = P(II)
  192 CONTINUE
      IP = IP-MR
      GO TO 195
  193 DO 194 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)-Q(I,JM2)
  194 CONTINUE
  195 CONTINUE
      CALL COSGEN (KR,1,FNUM2,0.5,TCOS)
      CALL COSGEN (LR,1,FNUM2,0.5,TCOS(KR+1))
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 196 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)+B(I)
  196 CONTINUE
      NLASTP = NLAST
  197 CONTINUE
      JSTEP = JR
      JR = I2R
      I2R = I2R/2
      IF (JR .EQ. 0) GO TO 210
      JSTART = 1+JR
      KR = KR-JR
      IF (NLAST+JR .GT. N) GO TO 198
      KR = KR-JR
      NLAST = NLAST+JR
      JSTOP = NLAST-JSTEP
      GO TO 199
  198 CONTINUE
      JSTOP = NLAST-JR
  199 CONTINUE
      LR = KR-JR
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      DO 209 J=JSTART,JSTOP,JSTEP
         JM2 = J-JR
         JP2 = J+JR
         IF (J .NE. JR) GO TO 201
         DO 200 I=1,MR
            B(I) = Q(I,J)+Q(I,JP2)
  200    CONTINUE
         GO TO 203
  201    CONTINUE
         DO 202 I=1,MR
            B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2)
  202    CONTINUE
  203    CONTINUE
         IF (JR .NE. 1) GO TO 205
         DO 204 I=1,MR
            Q(I,J) = 0.
  204    CONTINUE
         GO TO 207
  205    CONTINUE
         JM1 = J-I2R
         JP1 = J+I2R
         DO 206 I=1,MR
            Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  206    CONTINUE
  207    CONTINUE
         CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
         DO 208 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  208    CONTINUE
  209 CONTINUE
      NROD = 1
      IF (NLAST+I2R .LE. N) NROD = 0
      IF (NLASTP .NE. NLAST) GO TO 188
      GO TO 197
  210 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
*DECK TRI3
      SUBROUTINE TRI3 (M, A, B, C, K, Y1, Y2, Y3, TCOS, D, W1, W2, W3)
C***BEGIN PROLOGUE  TRI3
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (TRI3-S, CMPTR3-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve three linear systems whose common coefficient
C     matrix is a rational function in the matrix given by
C
C                  TRIDIAGONAL (...,A(I),B(I),C(I),...)
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890206  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  TRI3
      DIMENSION       A(*)       ,B(*)       ,C(*)       ,K(4)       ,
     1                TCOS(*)    ,Y1(*)      ,Y2(*)      ,Y3(*)      ,
     2                D(*)       ,W1(*)      ,W2(*)      ,W3(*)
      INTEGER K1P1, K2P1, K3P1, K4P1
C
C***FIRST EXECUTABLE STATEMENT  TRI3
      MM1 = M-1
      K1 = K(1)
      K2 = K(2)
      K3 = K(3)
      K4 = K(4)
      K1P1 = K1+1
      K2P1 = K2+1
      K3P1 = K3+1
      K4P1 = K4+1
      K2K3K4 = K2+K3+K4
      IF (K2K3K4 .EQ. 0) GO TO 101
      L1 = (K1+1)/(K2+1)
      L2 = (K1+1)/(K3+1)
      L3 = (K1+1)/(K4+1)
      LINT1 = 1
      LINT2 = 1
      LINT3 = 1
      KINT1 = K1
      KINT2 = KINT1+K2
      KINT3 = KINT2+K3
  101 CONTINUE
      DO 115 N=1,K1
         X = TCOS(N)
         IF (K2K3K4 .EQ. 0) GO TO 107
         IF (N .NE. L1) GO TO 103
         DO 102 I=1,M
            W1(I) = Y1(I)
  102    CONTINUE
  103    IF (N .NE. L2) GO TO 105
         DO 104 I=1,M
            W2(I) = Y2(I)
  104    CONTINUE
  105    IF (N .NE. L3) GO TO 107
         DO 106 I=1,M
            W3(I) = Y3(I)
  106    CONTINUE
  107    CONTINUE
         Z = 1./(B(1)-X)
         D(1) = C(1)*Z
         Y1(1) = Y1(1)*Z
         Y2(1) = Y2(1)*Z
         Y3(1) = Y3(1)*Z
         DO 108 I=2,M
            Z = 1./(B(I)-X-A(I)*D(I-1))
            D(I) = C(I)*Z
            Y1(I) = (Y1(I)-A(I)*Y1(I-1))*Z
            Y2(I) = (Y2(I)-A(I)*Y2(I-1))*Z
            Y3(I) = (Y3(I)-A(I)*Y3(I-1))*Z
  108    CONTINUE
         DO 109 IP=1,MM1
            I = M-IP
            Y1(I) = Y1(I)-D(I)*Y1(I+1)
            Y2(I) = Y2(I)-D(I)*Y2(I+1)
            Y3(I) = Y3(I)-D(I)*Y3(I+1)
  109    CONTINUE
         IF (K2K3K4 .EQ. 0) GO TO 115
         IF (N .NE. L1) GO TO 111
         I = LINT1+KINT1
         XX = X-TCOS(I)
         DO 110 I=1,M
            Y1(I) = XX*Y1(I)+W1(I)
  110    CONTINUE
         LINT1 = LINT1+1
         L1 = (LINT1*K1P1)/K2P1
  111    IF (N .NE. L2) GO TO 113
         I = LINT2+KINT2
         XX = X-TCOS(I)
         DO 112 I=1,M
            Y2(I) = XX*Y2(I)+W2(I)
  112    CONTINUE
         LINT2 = LINT2+1
         L2 = (LINT2*K1P1)/K3P1
  113    IF (N .NE. L3) GO TO 115
         I = LINT3+KINT3
         XX = X-TCOS(I)
         DO 114 I=1,M
            Y3(I) = XX*Y3(I)+W3(I)
  114    CONTINUE
         LINT3 = LINT3+1
         L3 = (LINT3*K1P1)/K4P1
  115 CONTINUE
      RETURN
      END
*DECK TRIX
      SUBROUTINE TRIX (IDEGBR, IDEGCR, M, A, B, C, Y, TCOS, D, W)
C***BEGIN PROLOGUE  TRIX
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (TRIX-S, CMPTRX-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve a system of linear equations where the
C     coefficient matrix is a rational function in the matrix given by
C     TRIDIAGONAL  ( . . . , A(I), B(I), C(I), . . . ).
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  TRIX
C
      DIMENSION       A(*)       ,B(*)       ,C(*)       ,Y(*)       ,
     1                TCOS(*)    ,D(*)       ,W(*)
      INTEGER KB, KC
C***FIRST EXECUTABLE STATEMENT  TRIX
      MM1 = M-1
      KB = IDEGBR+1
      KC = IDEGCR+1
      L = (IDEGBR+1)/(IDEGCR+1)
      LINT = 1
      DO 108 K=1,IDEGBR
         X = TCOS(K)
         IF (K .NE. L) GO TO 102
         I = IDEGBR+LINT
         XX = X-TCOS(I)
         DO 101 I=1,M
            W(I) = Y(I)
            Y(I) = XX*Y(I)
  101    CONTINUE
  102    CONTINUE
         Z = 1./(B(1)-X)
         D(1) = C(1)*Z
         Y(1) = Y(1)*Z
         DO 103 I=2,MM1
            Z = 1./(B(I)-X-A(I)*D(I-1))
            D(I) = C(I)*Z
            Y(I) = (Y(I)-A(I)*Y(I-1))*Z
  103    CONTINUE
         Z = B(M)-X-A(M)*D(MM1)
         IF (Z .NE. 0.) GO TO 104
         Y(M) = 0.
         GO TO 105
  104    Y(M) = (Y(M)-A(M)*Y(MM1))/Z
  105    CONTINUE
         DO 106 IP=1,MM1
            I = M-IP
            Y(I) = Y(I)-D(I)*Y(I+1)
  106    CONTINUE
         IF (K .NE. L) GO TO 108
         DO 107 I=1,M
            Y(I) = Y(I)+W(I)
  107    CONTINUE
         LINT = LINT+1
         L = (LINT*KB)/KC
  108 CONTINUE
      RETURN
      END
*DECK S1MERG
      SUBROUTINE S1MERG (TCOS, I1, M1, I2, M2, I3)
C***BEGIN PROLOGUE  S1MERG
C***SUBSIDIARY
C***PURPOSE  Merge two strings of ascending real numbers.
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (S1MERG-S, D1MERG-D, C1MERG-C, I1MERG-I)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   This subroutine merges two ascending strings of numbers in the
C   array TCOS.  The first string is of length M1 and starts at
C   TCOS(I1+1).  The second string is of length M2 and starts at
C   TCOS(I2+1).  The merged string goes into TCOS(I3+1).
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  SCOPY
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   901120  Modified to use IF-THEN-ELSE.  Previous spaghetti code did
C           not compile correctly with optimization on the IBM RS6000.
C           (RWC)
C   920130  Code name changed from MERGE to S1MERG.  (WRB)
C***END PROLOGUE  S1MERG
      INTEGER I1, I2, I3, M1, M2
      REAL TCOS(*)
C
      INTEGER J1, J2, J3
C
C***FIRST EXECUTABLE STATEMENT  S1MERG
      IF (M1.EQ.0 .AND. M2.EQ.0) RETURN
C
      IF (M1.EQ.0 .AND. M2.NE.0) THEN
         CALL SCOPY (M2, TCOS(I2+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      IF (M1.NE.0 .AND. M2.EQ.0) THEN
         CALL SCOPY (M1, TCOS(I1+1), 1, TCOS(I3+1), 1)
         RETURN
      ENDIF
C
      J1 = 1
      J2 = 1
      J3 = 1
C
   10 IF (TCOS(I1+J1) .LE. TCOS(I2+J2)) THEN
         TCOS(I3+J3) = TCOS(I1+J1)
         J1 = J1+1
         IF (J1 .GT. M1) THEN
            CALL SCOPY (M2-J2+1, TCOS(I2+J2), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ELSE
         TCOS(I3+J3) = TCOS(I2+J2)
         J2 = J2+1
         IF (J2 .GT. M2) THEN
            CALL SCOPY (M1-J1+1, TCOS(I1+J1), 1, TCOS(I3+J3+1), 1)
            RETURN
         ENDIF
      ENDIF
      J3 = J3+1
      GO TO 10
      END
*DECK SCOPY
      SUBROUTINE SCOPY (N, SX, INCX, SY, INCY)
C***BEGIN PROLOGUE  SCOPY
C***PURPOSE  Copy a vector.
C***LIBRARY   SLATEC (BLAS)
C***CATEGORY  D1A5
C***TYPE      SINGLE PRECISION (SCOPY-S, DCOPY-D, CCOPY-C, ICOPY-I)
C***KEYWORDS  BLAS, COPY, LINEAR ALGEBRA, VECTOR
C***AUTHOR  Lawson, C. L., (JPL)
C           Hanson, R. J., (SNLA)
C           Kincaid, D. R., (U. of Texas)
C           Krogh, F. T., (JPL)
C***DESCRIPTION
C
C                B L A S  Subprogram
C    Description of Parameters
C
C     --Input--
C        N  number of elements in input vector(s)
C       SX  single precision vector with N elements
C     INCX  storage spacing between elements of SX
C       SY  single precision vector with N elements
C     INCY  storage spacing between elements of SY
C
C     --Output--
C       SY  copy of vector SX (unchanged if N .LE. 0)
C
C     Copy single precision SX to single precision SY.
C     For I = 0 to N-1, copy  SX(LX+I*INCX) to SY(LY+I*INCY),
C     where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
C     defined in a similar way using INCY.
C
C***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
C                 Krogh, Basic linear algebra subprograms for Fortran
C                 usage, Algorithm No. 539, Transactions on Mathematical
C                 Software 5, 3 (September 1979), pp. 308-323.
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   791001  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920310  Corrected definition of LX in DESCRIPTION.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  SCOPY
      REAL SX(*), SY(*)
C***FIRST EXECUTABLE STATEMENT  SCOPY
      IF (N .LE. 0) RETURN
      IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
C
C     Code for unequal or nonpositive increments.
C
    5 IX = 1
      IY = 1
      IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
      IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
      DO 10 I = 1,N
        SY(IY) = SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C     Code for both increments equal to 1.
C
C     Clean-up loop so remaining vector length is a multiple of 7.
C
   20 M = MOD(N,7)
      IF (M .EQ. 0) GO TO 40
      DO 30 I = 1,M
        SY(I) = SX(I)
   30 CONTINUE
      IF (N .LT. 7) RETURN
   40 MP1 = M + 1
      DO 50 I = MP1,N,7
        SY(I) = SX(I)
        SY(I+1) = SX(I+1)
        SY(I+2) = SX(I+2)
        SY(I+3) = SX(I+3)
        SY(I+4) = SX(I+4)
        SY(I+5) = SX(I+5)
        SY(I+6) = SX(I+6)
   50 CONTINUE
      RETURN
C
C     Code for equal, positive, non-unit increments.
C
   60 NS = N*INCX
      DO 70 I = 1,NS,INCX
        SY(I) = SX(I)
   70 CONTINUE
      RETURN
      END
*DECK GENBUN
      SUBROUTINE GENBUN (NPEROD, N, MPEROD, M, A, B, C, IDIMY, Y,
     +   IERROR, W)
C***BEGIN PROLOGUE  GENBUN
C***PURPOSE  Solve by a cyclic reduction algorithm the linear system
C            of equations that results from a finite difference
C            approximation to certain 2-d elliptic PDE's on a centered
C            grid .
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B4B
C***TYPE      SINGLE PRECISION (GENBUN-S, CMGNBN-C)
C***KEYWORDS  ELLIPTIC, FISHPACK, PDE, TRIDIAGONAL
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine GENBUN solves the linear system of equations
C
C          A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
C
C          + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
C
C               for I = 1,2,...,M  and  J = 1,2,...,N.
C
C     The indices I+1 and I-1 are evaluated modulo M, i.e.,
C     X(0,J) = X(M,J) and X(M+1,J) = X(1,J), and X(I,0) may be equal to
C     0, X(I,2), or X(I,N) and X(I,N+1) may be equal to 0, X(I,N-1), or
C     X(I,1) depending on an input parameter.
C
C
C     * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C             * * * * * *   On Input    * * * * * *
C
C     NPEROD
C       Indicates the values that X(I,0) and X(I,N+1) are assumed to
C       have.
C
C       = 0  If X(I,0) = X(I,N) and X(I,N+1) = X(I,1).
C       = 1  If X(I,0) = X(I,N+1) = 0  .
C       = 2  If X(I,0) = 0 and X(I,N+1) = X(I,N-1).
C       = 3  If X(I,0) = X(I,2) and X(I,N+1) = X(I,N-1).
C       = 4  If X(I,0) = X(I,2) and X(I,N+1) = 0.
C
C     N
C       The number of unknowns in the J-direction.  N must be greater
C       than 2.
C
C     MPEROD
C       = 0 if A(1) and C(M) are not zero.
C       = 1 if A(1) = C(M) = 0.
C
C     M
C       The number of unknowns in the I-direction.  M must be greater
C       than 2.
C
C     A,B,C
C       One-dimensional arrays of length M that specify the
C       coefficients in the linear equations given above.  If MPEROD = 0
C       the array elements must not depend upon the index I, but must be
C       constant.  Specifically, the subroutine checks the following
C       condition
C
C             A(I) = C(1)
C             C(I) = C(1)
C             B(I) = B(1)
C
C       for I=1,2,...,M.
C
C     IDIMY
C       The row (or first) dimension of the two-dimensional array Y as
C       it appears in the program calling GENBUN.  This parameter is
C       used to specify the variable dimension of Y.  IDIMY must be at
C       least M.
C
C     Y
C       A two-dimensional array that specifies the values of the right
C       side of the linear system of equations given above.  Y must be
C       dimensioned at least M*N.
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 + (10 + INT(log2(N)))*M
C       locations.  The actual number of locations used is computed by
C       GENBUN and is returned in location W(1).
C
C
C             * * * * * *   On Output     * * * * * *
C
C     Y
C       Contains the solution X.
C
C     IERROR
C       An error flag that indicates invalid input parameters.  Except
C       for number zero, a solution is not attempted.
C
C       = 0  No error.
C       = 1  M .LE. 2
C       = 2  N .LE. 2
C       = 3  IDIMY .LT. M
C       = 4  NPEROD .LT. 0 or NPEROD .GT. 4
C       = 5  MPEROD .LT. 0 or MPEROD .GT. 1
C       = 6  A(I) .NE. C(1) or C(I) .NE. C(1) or B(I) .NE. B(1) for
C            some I=1,2,...,M.
C       = 7  A(1) .NE. 0 or C(M) .NE. 0 and MPEROD = 1
C
C     W
C       W(1) contains the required length of W.
C
C *Long Description:
C
C     * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   A(M),B(M),C(M),Y(IDIMY,N),W(see parameter list)
C     Arguments
C
C     Latest         June 1, 1976
C     Revision
C
C     Subprograms    GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE,TRIX,TRI3,
C     Required       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     Roland Sweet
C
C     Language       FORTRAN
C
C     History        Standardized April 1, 1973
C                    Revised August 20,1973
C                    Revised January 1, 1976
C
C     Algorithm      The linear system is solved by a cyclic reduction
C                    algorithm described in the reference.
C
C     Space          4944(decimal) = 11520(octal) locations on the NCAR
C     Required       Control Data 7600.
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine GENBUN is roughly proportional
C                    to M*N*log2(N), but also depends on the input
C                    parameter NPEROD.  Some typical values are listed
C                    in the table below.  More comprehensive timing
C                    charts may be found in the reference.
C                       To measure the accuracy of the algorithm a
C                    uniform random number generator was used to create
C                    a solution array X for the system given in the
C                    'PURPOSE' with
C
C                       A(I) = C(I) = -0.5*B(I) = 1,       I=1,2,...,M
C
C                    and, when MPEROD = 1
C
C                       A(1) = C(M) = 0
C                       A(M) = C(1) = 2.
C
C                    The solution X was substituted into the given sys-
C                    tem and, using double precision, a right side Y was
C                    computed.  Using this array Y subroutine GENBUN was
C                    called to produce an approximate solution Z.  Then
C                    the relative error, defined as
C
C                       E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J)))
C
C                    where the two maxima are taken over all I=1,2,...,M
C                    and J=1,2,...,N, was computed.  The value of E is
C                    given in the table below for some typical values of
C                    M and N.
C
C
C                       M (=N)    MPEROD    NPEROD    T(MSECS)    E
C                       ------    ------    ------    --------  ------
C
C                         31        0         0          36     6.E-14
C                         31        1         1          21     4.E-13
C                         31        1         3          41     3.E-13
C                         32        0         0          29     9.E-14
C                         32        1         1          32     3.E-13
C                         32        1         3          48     1.E-13
C                         33        0         0          36     9.E-14
C                         33        1         1          30     4.E-13
C                         33        1         3          34     1.E-13
C                         63        0         0         150     1.E-13
C                         63        1         1          91     1.E-12
C                         63        1         3         173     2.E-13
C                         64        0         0         122     1.E-13
C                         64        1         1         128     1.E-12
C                         64        1         3         199     6.E-13
C                         65        0         0         143     2.E-13
C                         65        1         1         120     1.E-12
C                         65        1         3         138     4.E-13
C
C     Portability    American National Standards Institute Fortran.
C                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       COS
C     Resident
C     Routines
C
C     Reference      Sweet, R., 'A Cyclic Reduction Algorithm For
C                    Solving Block Tridiagonal Systems Of Arbitrary
C                    Dimensions,' SIAM J. on Numer. Anal.,
C                    14(Sept., 1977), PP. 706-720.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES  R. Sweet, A cyclic reduction algorithm for solving
C                 block tridiagonal systems of arbitrary dimensions,
C                 SIAM Journal on Numerical Analysis 14, (September
C                 1977), pp. 706-720.
C***ROUTINES CALLED  POISD2, POISN2, POISP2
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   861211  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  GENBUN
C
C
      DIMENSION       Y(IDIMY,*)
      DIMENSION       W(*)       ,B(*)       ,A(*)       ,C(*)
C***FIRST EXECUTABLE STATEMENT  GENBUN
      IERROR = 0
      IF (M .LE. 2) IERROR = 1
      IF (N .LE. 2) IERROR = 2
      IF (IDIMY .LT. M) IERROR = 3
      IF (NPEROD.LT.0 .OR. NPEROD.GT.4) IERROR = 4
      IF (MPEROD.LT.0 .OR. MPEROD.GT.1) IERROR = 5
      IF (MPEROD .EQ. 1) GO TO 102
      DO 101 I=2,M
         IF (A(I) .NE. C(1)) GO TO 103
         IF (C(I) .NE. C(1)) GO TO 103
         IF (B(I) .NE. B(1)) GO TO 103
  101 CONTINUE
      GO TO 104
  102 IF (A(1).NE.0. .OR. C(M).NE.0.) IERROR = 7
      GO TO 104
  103 IERROR = 6
  104 IF (IERROR .NE. 0) RETURN
      MP1 = M+1
      IWBA = MP1
      IWBB = IWBA+M
      IWBC = IWBB+M
      IWB2 = IWBC+M
      IWB3 = IWB2+M
      IWW1 = IWB3+M
      IWW2 = IWW1+M
      IWW3 = IWW2+M
      IWD = IWW3+M
      IWTCOS = IWD+M
      IWP = IWTCOS+4*N
      DO 106 I=1,M
         K = IWBA+I-1
         W(K) = -A(I)
         K = IWBC+I-1
         W(K) = -C(I)
         K = IWBB+I-1
         W(K) = 2.-B(I)
         DO 105 J=1,N
            Y(I,J) = -Y(I,J)
  105    CONTINUE
  106 CONTINUE
      MP = MPEROD+1
      NP = NPEROD+1
      GO TO (114,107),MP
  107 GO TO (108,109,110,111,123),NP
  108 CALL POISP2 (M,N,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
      GO TO 112
  109 CALL POISD2 (M,N,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWW1),
     1             W(IWD),W(IWTCOS),W(IWP))
      GO TO 112
  110 CALL POISN2 (M,N,1,2,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
      GO TO 112
  111 CALL POISN2 (M,N,1,1,W(IWBA),W(IWBB),W(IWBC),Y,IDIMY,W,W(IWB2),
     1             W(IWB3),W(IWW1),W(IWW2),W(IWW3),W(IWD),W(IWTCOS),
     2             W(IWP))
  112 IPSTOR = W(IWW1)
      IREV = 2
      IF (NPEROD .EQ. 4) GO TO 124
  113 GO TO (127,133),MP
  114 CONTINUE
C
C     REORDER UNKNOWNS WHEN MP =0
C
      MH = (M+1)/2
      MHM1 = MH-1
      MODD = 1
      IF (MH*2 .EQ. M) MODD = 2
      DO 119 J=1,N
         DO 115 I=1,MHM1
            MHPI = MH+I
            MHMI = MH-I
            W(I) = Y(MHMI,J)-Y(MHPI,J)
            W(MHPI) = Y(MHMI,J)+Y(MHPI,J)
  115    CONTINUE
         W(MH) = 2.*Y(MH,J)
         GO TO (117,116),MODD
  116    W(M) = 2.*Y(M,J)
  117    CONTINUE
         DO 118 I=1,M
            Y(I,J) = W(I)
  118    CONTINUE
  119 CONTINUE
      K = IWBC+MHM1-1
      I = IWBA+MHM1
      W(K) = 0.
      W(I) = 0.
      W(K+1) = 2.*W(K+1)
      GO TO (120,121),MODD
  120 CONTINUE
      K = IWBB+MHM1-1
      W(K) = W(K)-W(I-1)
      W(IWBC-1) = W(IWBC-1)+W(IWBB-1)
      GO TO 122
  121 W(IWBB-1) = W(K+1)
  122 CONTINUE
      GO TO 107
C
C     REVERSE COLUMNS WHEN NPEROD = 4.
C
  123 IREV = 1
      NBY2 = N/2
  124 DO 126 J=1,NBY2
         MSKIP = N+1-J
         DO 125 I=1,M
            A1 = Y(I,J)
            Y(I,J) = Y(I,MSKIP)
            Y(I,MSKIP) = A1
  125    CONTINUE
  126 CONTINUE
      GO TO (110,113),IREV
  127 CONTINUE
      DO 132 J=1,N
         DO 128 I=1,MHM1
            MHMI = MH-I
            MHPI = MH+I
            W(MHMI) = .5*(Y(MHPI,J)+Y(I,J))
            W(MHPI) = .5*(Y(MHPI,J)-Y(I,J))
  128    CONTINUE
         W(MH) = .5*Y(MH,J)
         GO TO (130,129),MODD
  129    W(M) = .5*Y(M,J)
  130    CONTINUE
         DO 131 I=1,M
            Y(I,J) = W(I)
  131    CONTINUE
  132 CONTINUE
  133 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR W ARRAY.
C
      W(1) = IPSTOR+IWP-1
      RETURN
      END
*DECK COSGEN
      SUBROUTINE COSGEN (N, IJUMP, FNUM, FDEN, A)
C***BEGIN PROLOGUE  COSGEN
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (COSGEN-S, CMPCSG-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     This subroutine computes required cosine values in ascending
C     order.  When IJUMP .GT. 1 the routine computes values
C
C        2*COS(J*PI/L) , J=1,2,...,L and J .NE. 0(MOD N/IJUMP+1)
C
C     where L = IJUMP*(N/IJUMP+1).
C
C
C     when IJUMP = 1 it computes
C
C            2*COS((J-FNUM)*PI/(N+FDEN)) ,  J=1, 2, ... ,N
C
C     where
C        FNUM = 0.5, FDEN = 0.0, for regular reduction values.
C        FNUM = 0.0, FDEN = 1.0, for B-R and C-R when ISTAG = 1
C        FNUM = 0.0, FDEN = 0.5, for B-R and C-R when ISTAG = 2
C        FNUM = 0.5, FDEN = 0.5, for B-R and C-R when ISTAG = 2
C                                in POISN2 only.
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  PIMACH
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  COSGEN
      DIMENSION       A(*)
C
C
C***FIRST EXECUTABLE STATEMENT  COSGEN
      PI = PIMACH(DUM)
      IF (N .EQ. 0) GO TO 105
      IF (IJUMP .EQ. 1) GO TO 103
      K3 = N/IJUMP+1
      K4 = K3-1
      PIBYN = PI/(N+IJUMP)
      DO 102 K=1,IJUMP
         K1 = (K-1)*K3
         K5 = (K-1)*K4
         DO 101 I=1,K4
            X = K1+I
            K2 = K5+I
            A(K2) = -2.*COS(X*PIBYN)
  101    CONTINUE
  102 CONTINUE
      GO TO 105
  103 CONTINUE
      NP1 = N+1
      Y = PI/(N+FDEN)
      DO 104 I=1,N
         X = NP1-I-FNUM
         A(I) = 2.*COS(X*Y)
  104 CONTINUE
  105 CONTINUE
      RETURN
      END
*DECK PIMACH
      FUNCTION PIMACH (DUM)
C***BEGIN PROLOGUE  PIMACH
C***SUBSIDIARY
C***PURPOSE  Subsidiary to HSTCSP, HSTSSP and HWSCSP
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (PIMACH-S)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     This subprogram supplies the value of the constant PI correct to
C     machine precision where
C
C     PI=3.1415926535897932384626433832795028841971693993751058209749446
C
C***SEE ALSO  HSTCSP, HSTSSP, HWSCSP
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  PIMACH
C
C***FIRST EXECUTABLE STATEMENT  PIMACH
      PIMACH = 3.14159265358979
      RETURN
      END
*DECK POISD2
      SUBROUTINE POISD2 (MR, NR, ISTAG, BA, BB, BC, Q, IDIMQ, B, W, D,
     +   TCOS, P)
C***BEGIN PROLOGUE  POISD2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (POISD2-S, CMPOSD-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve Poisson's equation for Dirichlet boundary
C     conditions.
C
C     ISTAG = 1 if the last diagonal block is the matrix A.
C     ISTAG = 2 if the last diagonal block is the matrix A+I.
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  COSGEN, S1MERG, TRIX
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   920130  Modified to use merge routine S1MERG rather than deleted
C           routine MERGE.  (WRB)
C***END PROLOGUE  POISD2
C
      DIMENSION       Q(IDIMQ,*) ,BA(*)      ,BB(*)      ,BC(*)      ,
     1                TCOS(*)    ,B(*)       ,D(*)       ,W(*)       ,
     2                P(*)
C***FIRST EXECUTABLE STATEMENT  POISD2
      M = MR
      N = NR
      JSH = 0
      FI = 1./ISTAG
      IP = -M
      IPSTOR = 0
      GO TO (101,102),ISTAG
  101 KR = 0
      IRREG = 1
      IF (N .GT. 1) GO TO 106
      TCOS(1) = 0.
      GO TO 103
  102 KR = 1
      JSTSAV = 1
      IRREG = 2
      IF (N .GT. 1) GO TO 106
      TCOS(1) = -1.
  103 DO 104 I=1,M
         B(I) = Q(I,1)
  104 CONTINUE
      CALL TRIX (1,0,M,BA,BB,BC,B,TCOS,D,W)
      DO 105 I=1,M
         Q(I,1) = B(I)
  105 CONTINUE
      GO TO 183
  106 LR = 0
      DO 107 I=1,M
         P(I) = 0.
  107 CONTINUE
      NUN = N
      JST = 1
      JSP = N
C
C     IRREG = 1 WHEN NO IRREGULARITIES HAVE OCCURRED, OTHERWISE IT IS 2.
C
  108 L = 2*JST
      NODD = 2-2*((NUN+1)/2)+NUN
C
C     NODD = 1 WHEN NUN IS ODD, OTHERWISE IT IS 2.
C
      GO TO (110,109),NODD
  109 JSP = JSP-L
      GO TO 111
  110 JSP = JSP-JST
      IF (IRREG .NE. 1) JSP = JSP-L
  111 CONTINUE
C
C     REGULAR REDUCTION
C
      CALL COSGEN (JST,1,0.5,0.0,TCOS)
      IF (L .GT. JSP) GO TO 118
      DO 117 J=L,JSP,L
         JM1 = J-JSH
         JP1 = J+JSH
         JM2 = J-JST
         JP2 = J+JST
         JM3 = JM2-JSH
         JP3 = JP2+JSH
         IF (JST .NE. 1) GO TO 113
         DO 112 I=1,M
            B(I) = 2.*Q(I,J)
            Q(I,J) = Q(I,JM2)+Q(I,JP2)
  112    CONTINUE
         GO TO 115
  113    DO 114 I=1,M
            T = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2)
            B(I) = T+Q(I,J)-Q(I,JM3)-Q(I,JP3)
            Q(I,J) = T
  114    CONTINUE
  115    CONTINUE
         CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W)
         DO 116 I=1,M
            Q(I,J) = Q(I,J)+B(I)
  116    CONTINUE
  117 CONTINUE
C
C     REDUCTION FOR LAST UNKNOWN
C
  118 GO TO (119,136),NODD
  119 GO TO (152,120),IRREG
C
C     ODD NUMBER OF UNKNOWNS
C
  120 JSP = JSP+L
      J = JSP
      JM1 = J-JSH
      JP1 = J+JSH
      JM2 = J-JST
      JP2 = J+JST
      JM3 = JM2-JSH
      GO TO (123,121),ISTAG
  121 CONTINUE
      IF (JST .NE. 1) GO TO 123
      DO 122 I=1,M
         B(I) = Q(I,J)
         Q(I,J) = 0.
  122 CONTINUE
      GO TO 130
  123 GO TO (124,126),NODDPR
  124 DO 125 I=1,M
         IP1 = IP+I
         B(I) = .5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+P(IP1)+Q(I,J)
  125 CONTINUE
      GO TO 128
  126 DO 127 I=1,M
         B(I) = .5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))+Q(I,JP2)-Q(I,JP1)+Q(I,J)
  127 CONTINUE
  128 DO 129 I=1,M
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  129 CONTINUE
  130 CALL TRIX (JST,0,M,BA,BB,BC,B,TCOS,D,W)
      IP = IP+M
      IPSTOR = MAX(IPSTOR,IP+M)
      DO 131 I=1,M
         IP1 = IP+I
         P(IP1) = Q(I,J)+B(I)
         B(I) = Q(I,JP2)+P(IP1)
  131 CONTINUE
      IF (LR .NE. 0) GO TO 133
      DO 132 I=1,JST
         KRPI = KR+I
         TCOS(KRPI) = TCOS(I)
  132 CONTINUE
      GO TO 134
  133 CONTINUE
      CALL COSGEN (LR,JSTSAV,0.,FI,TCOS(JST+1))
      CALL S1MERG (TCOS,0,JST,JST,LR,KR)
  134 CONTINUE
      CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL TRIX (KR,KR,M,BA,BB,BC,B,TCOS,D,W)
      DO 135 I=1,M
         IP1 = IP+I
         Q(I,J) = Q(I,JM2)+B(I)+P(IP1)
  135 CONTINUE
      LR = KR
      KR = KR+L
      GO TO 152
C
C     EVEN NUMBER OF UNKNOWNS
C
  136 JSP = JSP+L
      J = JSP
      JM1 = J-JSH
      JP1 = J+JSH
      JM2 = J-JST
      JP2 = J+JST
      JM3 = JM2-JSH
      GO TO (137,138),IRREG
  137 CONTINUE
      JSTSAV = JST
      IDEG = JST
      KR = L
      GO TO 139
  138 CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
      IDEG = KR
      KR = KR+JST
  139 IF (JST .NE. 1) GO TO 141
      IRREG = 2
      DO 140 I=1,M
         B(I) = Q(I,J)
         Q(I,J) = Q(I,JM2)
  140 CONTINUE
      GO TO 150
  141 DO 142 I=1,M
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  142 CONTINUE
      GO TO (143,145),IRREG
  143 DO 144 I=1,M
         Q(I,J) = Q(I,JM2)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  144 CONTINUE
      IRREG = 2
      GO TO 150
  145 CONTINUE
      GO TO (146,148),NODDPR
  146 DO 147 I=1,M
         IP1 = IP+I
         Q(I,J) = Q(I,JM2)+P(IP1)
  147 CONTINUE
      IP = IP-M
      GO TO 150
  148 DO 149 I=1,M
         Q(I,J) = Q(I,JM2)+Q(I,J)-Q(I,JM1)
  149 CONTINUE
  150 CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W)
      DO 151 I=1,M
         Q(I,J) = Q(I,J)+B(I)
  151 CONTINUE
  152 NUN = NUN/2
      NODDPR = NODD
      JSH = JST
      JST = 2*JST
      IF (NUN .GE. 2) GO TO 108
C
C     START SOLUTION.
C
      J = JSP
      DO 153 I=1,M
         B(I) = Q(I,J)
  153 CONTINUE
      GO TO (154,155),IRREG
  154 CONTINUE
      CALL COSGEN (JST,1,0.5,0.0,TCOS)
      IDEG = JST
      GO TO 156
  155 KR = LR+JST
      CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
      CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
      IDEG = KR
  156 CONTINUE
      CALL TRIX (IDEG,LR,M,BA,BB,BC,B,TCOS,D,W)
      JM1 = J-JSH
      JP1 = J+JSH
      GO TO (157,159),IRREG
  157 DO 158 I=1,M
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
  158 CONTINUE
      GO TO 164
  159 GO TO (160,162),NODDPR
  160 DO 161 I=1,M
         IP1 = IP+I
         Q(I,J) = P(IP1)+B(I)
  161 CONTINUE
      IP = IP-M
      GO TO 164
  162 DO 163 I=1,M
         Q(I,J) = Q(I,J)-Q(I,JM1)+B(I)
  163 CONTINUE
  164 CONTINUE
C
C     START BACK SUBSTITUTION.
C
      JST = JST/2
      JSH = JST/2
      NUN = 2*NUN
      IF (NUN .GT. N) GO TO 183
      DO 182 J=JST,N,L
         JM1 = J-JSH
         JP1 = J+JSH
         JM2 = J-JST
         JP2 = J+JST
         IF (J .GT. JST) GO TO 166
         DO 165 I=1,M
            B(I) = Q(I,J)+Q(I,JP2)
  165    CONTINUE
         GO TO 170
  166    IF (JP2 .LE. N) GO TO 168
         DO 167 I=1,M
            B(I) = Q(I,J)+Q(I,JM2)
  167    CONTINUE
         IF (JST .LT. JSTSAV) IRREG = 1
         GO TO (170,171),IRREG
  168    DO 169 I=1,M
            B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2)
  169    CONTINUE
  170    CONTINUE
         CALL COSGEN (JST,1,0.5,0.0,TCOS)
         IDEG = JST
         JDEG = 0
         GO TO 172
  171    IF (J+L .GT. N) LR = LR-JST
         KR = JST+LR
         CALL COSGEN (KR,JSTSAV,0.0,FI,TCOS)
         CALL COSGEN (LR,JSTSAV,0.0,FI,TCOS(KR+1))
         IDEG = KR
         JDEG = LR
  172    CONTINUE
         CALL TRIX (IDEG,JDEG,M,BA,BB,BC,B,TCOS,D,W)
         IF (JST .GT. 1) GO TO 174
         DO 173 I=1,M
            Q(I,J) = B(I)
  173    CONTINUE
         GO TO 182
  174    IF (JP2 .GT. N) GO TO 177
  175    DO 176 I=1,M
            Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
  176    CONTINUE
         GO TO 182
  177    GO TO (175,178),IRREG
  178    IF (J+JSH .GT. N) GO TO 180
         DO 179 I=1,M
            IP1 = IP+I
            Q(I,J) = B(I)+P(IP1)
  179    CONTINUE
         IP = IP-M
         GO TO 182
  180    DO 181 I=1,M
            Q(I,J) = B(I)+Q(I,J)-Q(I,JM1)
  181    CONTINUE
  182 CONTINUE
      L = L/2
      GO TO 164
  183 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
*DECK POISN2
      SUBROUTINE POISN2 (M, N, ISTAG, MIXBND, A, BB, C, Q, IDIMQ, B, B2,
     +   B3, W, W2, W3, D, TCOS, P)
C***BEGIN PROLOGUE  POISN2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (POISN2-S, CMPOSN-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve Poisson's equation with Neumann boundary
C     conditions.
C
C     ISTAG = 1 if the last diagonal block is A.
C     ISTAG = 2 if the last diagonal block is A-I.
C     MIXBND = 1 if have Neumann boundary conditions at both boundaries.
C     MIXBND = 2 if have Neumann boundary conditions at bottom and
C     Dirichlet condition at top.  (for this case, must have ISTAG = 1.)
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  COSGEN, S1MERG, TRI3, TRIX
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C   920130  Modified to use merge routine S1MERG rather than deleted
C           routine MERGE.  (WRB)
C***END PROLOGUE  POISN2
C
      DIMENSION       A(*)       ,BB(*)      ,C(*)       ,Q(IDIMQ,*) ,
     1                B(*)       ,B2(*)      ,B3(*)      ,W(*)       ,
     2                W2(*)      ,W3(*)      ,D(*)       ,TCOS(*)    ,
     3                K(4)       ,P(*)
      EQUIVALENCE     (K(1),K1)  ,(K(2),K2)  ,(K(3),K3)  ,(K(4),K4)
C***FIRST EXECUTABLE STATEMENT  POISN2
      FISTAG = 3-ISTAG
      FNUM = 1./ISTAG
      FDEN = 0.5*(ISTAG-1)
      MR = M
      IP = -MR
      IPSTOR = 0
      I2R = 1
      JR = 2
      NR = N
      NLAST = N
      KR = 1
      LR = 0
      GO TO (101,103),ISTAG
  101 CONTINUE
      DO 102 I=1,MR
         Q(I,N) = .5*Q(I,N)
  102 CONTINUE
      GO TO (103,104),MIXBND
  103 IF (N .LE. 3) GO TO 155
  104 CONTINUE
      JR = 2*I2R
      NROD = 1
      IF ((NR/2)*2 .EQ. NR) NROD = 0
      GO TO (105,106),MIXBND
  105 JSTART = 1
      GO TO 107
  106 JSTART = JR
      NROD = 1-NROD
  107 CONTINUE
      JSTOP = NLAST-JR
      IF (NROD .EQ. 0) JSTOP = JSTOP-I2R
      CALL COSGEN (I2R,1,0.5,0.0,TCOS)
      I2RBY2 = I2R/2
      IF (JSTOP .GE. JSTART) GO TO 108
      J = JR
      GO TO 116
  108 CONTINUE
C
C     REGULAR REDUCTION.
C
      DO 115 J=JSTART,JSTOP,JR
         JP1 = J+I2RBY2
         JP2 = J+I2R
         JP3 = JP2+I2RBY2
         JM1 = J-I2RBY2
         JM2 = J-I2R
         JM3 = JM2-I2RBY2
         IF (J .NE. 1) GO TO 109
         JM1 = JP1
         JM2 = JP2
         JM3 = JP3
  109    CONTINUE
         IF (I2R .NE. 1) GO TO 111
         IF (J .EQ. 1) JM2 = JP2
         DO 110 I=1,MR
            B(I) = 2.*Q(I,J)
            Q(I,J) = Q(I,JM2)+Q(I,JP2)
  110    CONTINUE
         GO TO 113
  111    CONTINUE
         DO 112 I=1,MR
            FI = Q(I,J)
            Q(I,J) = Q(I,J)-Q(I,JM1)-Q(I,JP1)+Q(I,JM2)+Q(I,JP2)
            B(I) = FI+Q(I,J)-Q(I,JM3)-Q(I,JP3)
  112    CONTINUE
  113    CONTINUE
         CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
         DO 114 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  114    CONTINUE
C
C     END OF REDUCTION FOR REGULAR UNKNOWNS.
C
  115 CONTINUE
C
C     BEGIN SPECIAL REDUCTION FOR LAST UNKNOWN.
C
      J = JSTOP+JR
  116 NLAST = J
      JM1 = J-I2RBY2
      JM2 = J-I2R
      JM3 = JM2-I2RBY2
      IF (NROD .EQ. 0) GO TO 128
C
C     ODD NUMBER OF UNKNOWNS
C
      IF (I2R .NE. 1) GO TO 118
      DO 117 I=1,MR
         B(I) = FISTAG*Q(I,J)
         Q(I,J) = Q(I,JM2)
  117 CONTINUE
      GO TO 126
  118 DO 119 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  119 CONTINUE
      IF (NRODPR .NE. 0) GO TO 121
      DO 120 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)
  120 CONTINUE
      IP = IP-MR
      GO TO 123
  121 CONTINUE
      DO 122 I=1,MR
         Q(I,J) = Q(I,J)-Q(I,JM1)+Q(I,JM2)
  122 CONTINUE
  123 IF (LR .EQ. 0) GO TO 124
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1))
      GO TO 126
  124 CONTINUE
      DO 125 I=1,MR
         B(I) = FISTAG*B(I)
  125 CONTINUE
  126 CONTINUE
      CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 127 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
  127 CONTINUE
      KR = KR+I2R
      GO TO 151
  128 CONTINUE
C
C     EVEN NUMBER OF UNKNOWNS
C
      JP1 = J+I2RBY2
      JP2 = J+I2R
      IF (I2R .NE. 1) GO TO 135
      DO 129 I=1,MR
         B(I) = Q(I,J)
  129 CONTINUE
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      IP = 0
      IPSTOR = MR
      GO TO (133,130),ISTAG
  130 DO 131 I=1,MR
         P(I) = B(I)
         B(I) = B(I)+Q(I,N)
  131 CONTINUE
      TCOS(1) = 1.
      TCOS(2) = 0.
      CALL TRIX (1,1,MR,A,BB,C,B,TCOS,D,W)
      DO 132 I=1,MR
         Q(I,J) = Q(I,JM2)+P(I)+B(I)
  132 CONTINUE
      GO TO 150
  133 CONTINUE
      DO 134 I=1,MR
         P(I) = B(I)
         Q(I,J) = Q(I,JM2)+2.*Q(I,JP2)+3.*B(I)
  134 CONTINUE
      GO TO 150
  135 CONTINUE
      DO 136 I=1,MR
         B(I) = Q(I,J)+.5*(Q(I,JM2)-Q(I,JM1)-Q(I,JM3))
  136 CONTINUE
      IF (NRODPR .NE. 0) GO TO 138
      DO 137 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  137 CONTINUE
      GO TO 140
  138 CONTINUE
      DO 139 I=1,MR
         B(I) = B(I)+Q(I,JP2)-Q(I,JP1)
  139 CONTINUE
  140 CONTINUE
      CALL TRIX (I2R,0,MR,A,BB,C,B,TCOS,D,W)
      IP = IP+MR
      IPSTOR = MAX(IPSTOR,IP+MR)
      DO 141 I=1,MR
         II = IP+I
         P(II) = B(I)+.5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         B(I) = P(II)+Q(I,JP2)
  141 CONTINUE
      IF (LR .EQ. 0) GO TO 142
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(I2R+1))
      CALL S1MERG (TCOS,0,I2R,I2R,LR,KR)
      GO TO 144
  142 DO 143 I=1,I2R
         II = KR+I
         TCOS(II) = TCOS(I)
  143 CONTINUE
  144 CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      IF (LR .NE. 0) GO TO 145
      GO TO (146,145),ISTAG
  145 CONTINUE
      CALL TRIX (KR,KR,MR,A,BB,C,B,TCOS,D,W)
      GO TO 148
  146 CONTINUE
      DO 147 I=1,MR
         B(I) = FISTAG*B(I)
  147 CONTINUE
  148 CONTINUE
      DO 149 I=1,MR
         II = IP+I
         Q(I,J) = Q(I,JM2)+P(II)+B(I)
  149 CONTINUE
  150 CONTINUE
      LR = KR
      KR = KR+JR
  151 CONTINUE
      GO TO (152,153),MIXBND
  152 NR = (NLAST-1)/JR+1
      IF (NR .LE. 3) GO TO 155
      GO TO 154
  153 NR = NLAST/JR
      IF (NR .LE. 1) GO TO 192
  154 I2R = JR
      NRODPR = NROD
      GO TO 104
  155 CONTINUE
C
C      BEGIN SOLUTION
C
      J = 1+JR
      JM1 = J-I2R
      JP1 = J+I2R
      JM2 = NLAST-I2R
      IF (NR .EQ. 2) GO TO 184
      IF (LR .NE. 0) GO TO 170
      IF (N .NE. 3) GO TO 161
C
C     CASE N = 3.
C
      GO TO (156,168),ISTAG
  156 CONTINUE
      DO 157 I=1,MR
         B(I) = Q(I,2)
  157 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 158 I=1,MR
         Q(I,2) = B(I)
         B(I) = 4.*B(I)+Q(I,1)+2.*Q(I,3)
  158 CONTINUE
      TCOS(1) = -2.
      TCOS(2) = 2.
      I1 = 2
      I2 = 0
      CALL TRIX (I1,I2,MR,A,BB,C,B,TCOS,D,W)
      DO 159 I=1,MR
         Q(I,2) = Q(I,2)+B(I)
         B(I) = Q(I,1)+2.*Q(I,2)
  159 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 160 I=1,MR
         Q(I,1) = B(I)
  160 CONTINUE
      JR = 1
      I2R = 0
      GO TO 194
C
C     CASE N = 2**P+1
C
  161 CONTINUE
      GO TO (162,170),ISTAG
  162 CONTINUE
      DO 163 I=1,MR
         B(I) = Q(I,J)+.5*Q(I,1)-Q(I,JM1)+Q(I,NLAST)-Q(I,JM2)
  163 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      DO 164 I=1,MR
         Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))+B(I)
         B(I) = Q(I,1)+2.*Q(I,NLAST)+4.*Q(I,J)
  164 CONTINUE
      JR2 = 2*JR
      CALL COSGEN (JR,1,0.0,0.0,TCOS)
      DO 165 I=1,JR
         I1 = JR+I
         I2 = JR+1-I
         TCOS(I1) = -TCOS(I2)
  165 CONTINUE
      CALL TRIX (JR2,0,MR,A,BB,C,B,TCOS,D,W)
      DO 166 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
         B(I) = Q(I,1)+2.*Q(I,J)
  166 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      DO 167 I=1,MR
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)+B(I)
  167 CONTINUE
      GO TO 194
C
C     CASE OF GENERAL N WITH NR = 3 .
C
  168 DO 169 I=1,MR
         B(I) = Q(I,2)
         Q(I,2) = 0.
         B2(I) = Q(I,3)
         B3(I) = Q(I,1)
  169 CONTINUE
      JR = 1
      I2R = 0
      J = 2
      GO TO 177
  170 CONTINUE
      DO 171 I=1,MR
         B(I) = .5*Q(I,1)-Q(I,JM1)+Q(I,J)
  171 CONTINUE
      IF (NROD .NE. 0) GO TO 173
      DO 172 I=1,MR
         II = IP+I
         B(I) = B(I)+P(II)
  172 CONTINUE
      GO TO 175
  173 DO 174 I=1,MR
         B(I) = B(I)+Q(I,NLAST)-Q(I,JM2)
  174 CONTINUE
  175 CONTINUE
      DO 176 I=1,MR
         T = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
         Q(I,J) = T
         B2(I) = Q(I,NLAST)+T
         B3(I) = Q(I,1)+2.*T
  176 CONTINUE
  177 CONTINUE
      K1 = KR+2*JR-1
      K2 = KR+JR
      TCOS(K1+1) = -2.
      K4 = K1+3-ISTAG
      CALL COSGEN (K2+ISTAG-2,1,0.0,FNUM,TCOS(K4))
      K4 = K1+K2+1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K4))
      CALL S1MERG (TCOS,K1,K2,K1+K2,JR-1,0)
      K3 = K1+K2+LR
      CALL COSGEN (JR,1,0.5,0.0,TCOS(K3+1))
      K4 = K3+JR+1
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K4))
      CALL S1MERG (TCOS,K3,JR,K3+JR,KR,K1)
      IF (LR .EQ. 0) GO TO 178
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(K4))
      CALL S1MERG (TCOS,K3,JR,K3+JR,LR,K3-LR)
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K4))
  178 K3 = KR
      K4 = KR
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 179 I=1,MR
         B(I) = B(I)+B2(I)+B3(I)
  179 CONTINUE
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 180 I=1,MR
         Q(I,J) = Q(I,J)+B(I)
         B(I) = Q(I,1)+2.*Q(I,J)
  180 CONTINUE
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
      IF (JR .NE. 1) GO TO 182
      DO 181 I=1,MR
         Q(I,1) = B(I)
  181 CONTINUE
      GO TO 194
  182 CONTINUE
      DO 183 I=1,MR
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)+B(I)
  183 CONTINUE
      GO TO 194
  184 CONTINUE
      IF (N .NE. 2) GO TO 188
C
C     CASE  N = 2
C
      DO 185 I=1,MR
         B(I) = Q(I,1)
  185 CONTINUE
      TCOS(1) = 0.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 186 I=1,MR
         Q(I,1) = B(I)
         B(I) = 2.*(Q(I,2)+B(I))*FISTAG
  186 CONTINUE
      TCOS(1) = -FISTAG
      TCOS(2) = 2.
      CALL TRIX (2,0,MR,A,BB,C,B,TCOS,D,W)
      DO 187 I=1,MR
         Q(I,1) = Q(I,1)+B(I)
  187 CONTINUE
      JR = 1
      I2R = 0
      GO TO 194
  188 CONTINUE
C
C     CASE OF GENERAL N AND NR = 2 .
C
      DO 189 I=1,MR
         II = IP+I
         B3(I) = 0.
         B(I) = Q(I,1)+2.*P(II)
         Q(I,1) = .5*Q(I,1)-Q(I,JM1)
         B2(I) = 2.*(Q(I,1)+Q(I,NLAST))
  189 CONTINUE
      K1 = KR+JR-1
      TCOS(K1+1) = -2.
      K4 = K1+3-ISTAG
      CALL COSGEN (KR+ISTAG-2,1,0.0,FNUM,TCOS(K4))
      K4 = K1+KR+1
      CALL COSGEN (JR-1,1,0.0,1.0,TCOS(K4))
      CALL S1MERG (TCOS,K1,KR,K1+KR,JR-1,0)
      CALL COSGEN (KR,1,0.5,FDEN,TCOS(K1+1))
      K2 = KR
      K4 = K1+K2+1
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(K4))
      K3 = LR
      K4 = 0
      CALL TRI3 (MR,A,BB,C,K,B,B2,B3,TCOS,D,W,W2,W3)
      DO 190 I=1,MR
         B(I) = B(I)+B2(I)
  190 CONTINUE
      TCOS(1) = 2.
      CALL TRIX (1,0,MR,A,BB,C,B,TCOS,D,W)
      DO 191 I=1,MR
         Q(I,1) = Q(I,1)+B(I)
  191 CONTINUE
      GO TO 194
  192 DO 193 I=1,MR
         B(I) = Q(I,NLAST)
  193 CONTINUE
      GO TO 196
  194 CONTINUE
C
C     START BACK SUBSTITUTION.
C
      J = NLAST-JR
      DO 195 I=1,MR
         B(I) = Q(I,NLAST)+Q(I,J)
  195 CONTINUE
  196 JM2 = NLAST-I2R
      IF (JR .NE. 1) GO TO 198
      DO 197 I=1,MR
         Q(I,NLAST) = 0.
  197 CONTINUE
      GO TO 202
  198 CONTINUE
      IF (NROD .NE. 0) GO TO 200
      DO 199 I=1,MR
         II = IP+I
         Q(I,NLAST) = P(II)
  199 CONTINUE
      IP = IP-MR
      GO TO 202
  200 DO 201 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)-Q(I,JM2)
  201 CONTINUE
  202 CONTINUE
      CALL COSGEN (KR,1,0.5,FDEN,TCOS)
      CALL COSGEN (LR,1,0.5,FDEN,TCOS(KR+1))
      IF (LR .NE. 0) GO TO 204
      DO 203 I=1,MR
         B(I) = FISTAG*B(I)
  203 CONTINUE
  204 CONTINUE
      CALL TRIX (KR,LR,MR,A,BB,C,B,TCOS,D,W)
      DO 205 I=1,MR
         Q(I,NLAST) = Q(I,NLAST)+B(I)
  205 CONTINUE
      NLASTP = NLAST
  206 CONTINUE
      JSTEP = JR
      JR = I2R
      I2R = I2R/2
      IF (JR .EQ. 0) GO TO 222
      GO TO (207,208),MIXBND
  207 JSTART = 1+JR
      GO TO 209
  208 JSTART = JR
  209 CONTINUE
      KR = KR-JR
      IF (NLAST+JR .GT. N) GO TO 210
      KR = KR-JR
      NLAST = NLAST+JR
      JSTOP = NLAST-JSTEP
      GO TO 211
  210 CONTINUE
      JSTOP = NLAST-JR
  211 CONTINUE
      LR = KR-JR
      CALL COSGEN (JR,1,0.5,0.0,TCOS)
      DO 221 J=JSTART,JSTOP,JSTEP
         JM2 = J-JR
         JP2 = J+JR
         IF (J .NE. JR) GO TO 213
         DO 212 I=1,MR
            B(I) = Q(I,J)+Q(I,JP2)
  212    CONTINUE
         GO TO 215
  213    CONTINUE
         DO 214 I=1,MR
            B(I) = Q(I,J)+Q(I,JM2)+Q(I,JP2)
  214    CONTINUE
  215    CONTINUE
         IF (JR .NE. 1) GO TO 217
         DO 216 I=1,MR
            Q(I,J) = 0.
  216    CONTINUE
         GO TO 219
  217    CONTINUE
         JM1 = J-I2R
         JP1 = J+I2R
         DO 218 I=1,MR
            Q(I,J) = .5*(Q(I,J)-Q(I,JM1)-Q(I,JP1))
  218    CONTINUE
  219    CONTINUE
         CALL TRIX (JR,0,MR,A,BB,C,B,TCOS,D,W)
         DO 220 I=1,MR
            Q(I,J) = Q(I,J)+B(I)
  220    CONTINUE
  221 CONTINUE
      NROD = 1
      IF (NLAST+I2R .LE. N) NROD = 0
      IF (NLASTP .NE. NLAST) GO TO 194
      GO TO 206
  222 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
*DECK POISP2
      SUBROUTINE POISP2 (M, N, A, BB, C, Q, IDIMQ, B, B2, B3, W, W2, W3,
     +   D, TCOS, P)
C***BEGIN PROLOGUE  POISP2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to GENBUN
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (POISP2-S, CMPOSP-C)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     Subroutine to solve Poisson equation with periodic boundary
C     conditions.
C
C***SEE ALSO  GENBUN
C***ROUTINES CALLED  POISD2, POISN2
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900402  Added TYPE section.  (WRB)
C***END PROLOGUE  POISP2
C
      DIMENSION       A(*)       ,BB(*)      ,C(*)       ,Q(IDIMQ,*) ,
     1                B(*)       ,B2(*)      ,B3(*)      ,W(*)       ,
     2                W2(*)      ,W3(*)      ,D(*)       ,TCOS(*)    ,
     3                P(*)
C***FIRST EXECUTABLE STATEMENT  POISP2
      MR = M
      NR = (N+1)/2
      NRM1 = NR-1
      IF (2*NR .NE. N) GO TO 107
C
C     EVEN NUMBER OF UNKNOWNS
C
      DO 102 J=1,NRM1
         NRMJ = NR-J
         NRPJ = NR+J
         DO 101 I=1,MR
            S = Q(I,NRMJ)-Q(I,NRPJ)
            T = Q(I,NRMJ)+Q(I,NRPJ)
            Q(I,NRMJ) = S
            Q(I,NRPJ) = T
  101    CONTINUE
  102 CONTINUE
      DO 103 I=1,MR
         Q(I,NR) = 2.*Q(I,NR)
         Q(I,N) = 2.*Q(I,N)
  103 CONTINUE
      CALL POISD2 (MR,NRM1,1,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P)
      IPSTOR = W(1)
      CALL POISN2 (MR,NR+1,1,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D,
     1             TCOS,P)
      IPSTOR = MAX(IPSTOR,INT(W(1)))
      DO 105 J=1,NRM1
         NRMJ = NR-J
         NRPJ = NR+J
         DO 104 I=1,MR
            S = .5*(Q(I,NRPJ)+Q(I,NRMJ))
            T = .5*(Q(I,NRPJ)-Q(I,NRMJ))
            Q(I,NRMJ) = S
            Q(I,NRPJ) = T
  104    CONTINUE
  105 CONTINUE
      DO 106 I=1,MR
         Q(I,NR) = .5*Q(I,NR)
         Q(I,N) = .5*Q(I,N)
  106 CONTINUE
      GO TO 118
  107 CONTINUE
C
C     ODD  NUMBER OF UNKNOWNS
C
      DO 109 J=1,NRM1
         NRPJ = N+1-J
         DO 108 I=1,MR
            S = Q(I,J)-Q(I,NRPJ)
            T = Q(I,J)+Q(I,NRPJ)
            Q(I,J) = S
            Q(I,NRPJ) = T
  108    CONTINUE
  109 CONTINUE
      DO 110 I=1,MR
         Q(I,NR) = 2.*Q(I,NR)
  110 CONTINUE
      LH = NRM1/2
      DO 112 J=1,LH
         NRMJ = NR-J
         DO 111 I=1,MR
            S = Q(I,J)
            Q(I,J) = Q(I,NRMJ)
            Q(I,NRMJ) = S
  111    CONTINUE
  112 CONTINUE
      CALL POISD2 (MR,NRM1,2,A,BB,C,Q,IDIMQ,B,W,D,TCOS,P)
      IPSTOR = W(1)
      CALL POISN2 (MR,NR,2,1,A,BB,C,Q(1,NR),IDIMQ,B,B2,B3,W,W2,W3,D,
     1             TCOS,P)
      IPSTOR = MAX(IPSTOR,INT(W(1)))
      DO 114 J=1,NRM1
         NRPJ = NR+J
         DO 113 I=1,MR
            S = .5*(Q(I,NRPJ)+Q(I,J))
            T = .5*(Q(I,NRPJ)-Q(I,J))
            Q(I,NRPJ) = T
            Q(I,J) = S
  113    CONTINUE
  114 CONTINUE
      DO 115 I=1,MR
         Q(I,NR) = .5*Q(I,NR)
  115 CONTINUE
      DO 117 J=1,LH
         NRMJ = NR-J
         DO 116 I=1,MR
            S = Q(I,J)
            Q(I,J) = Q(I,NRMJ)
            Q(I,NRMJ) = S
  116    CONTINUE
  117 CONTINUE
  118 CONTINUE
C
C     RETURN STORAGE REQUIREMENTS FOR P VECTORS.
C
      W(1) = IPSTOR
      RETURN
      END
