# hw3crt.f

          SUBROUTINE HW3CRT (XS, XF, L, LBDCND, BDXS, BDXF, YS, YF, M,
         +   MBDCND, BDYS, BDYF, ZS, ZF, N, NBDCND, BDZS, BDZF, ELMBDA,
         +   LDIMF, MDIMF, F, PERTRB, IERROR, W)
    C***BEGIN PROLOGUE  HW3CRT
    C***PURPOSE  Solve the standard seven-point finite difference
    C            approximation to the Helmholtz equation in Cartesian
    C            coordinates.
    C***LIBRARY   SLATEC (FISHPACK)
    C***CATEGORY  I2B1A1A
    C***TYPE      SINGLE PRECISION (HW3CRT-S)
    C***KEYWORDS  CARTESIAN, ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
    C***AUTHOR  Adams, J., (NCAR)
    C           Swarztrauber, P. N., (NCAR)
    C           Sweet, R., (NCAR)
    C***DESCRIPTION
    C
    C     Subroutine HW3CRT solves the standard seven-point finite
    C     difference approximation to the Helmholtz equation in Cartesian
    C     coordinates:
    C
    C         (d/dX)(dU/dX) + (d/dY)(dU/dY) + (d/dZ)(dU/dZ)
    C
    C                    + LAMBDA*U = F(X,Y,Z) .
    C
    C    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
    C
    C
    C    * * * * * * * *    Parameter Description     * * * * * * * * * *
    C
    C
    C            * * * * * *   On Input    * * * * * *
    C
    C     XS,XF
    C        The range of X, i.e. XS .LE. X .LE. XF .
    C        XS must be less than XF.
    C
    C     L
    C        The number of panels into which the interval (XS,XF) is
    C        subdivided.  Hence, there will be L+1 grid points in the
    C        X-direction given by X(I) = XS+(I-1)DX for I=1,2,...,L+1,
    C        where DX = (XF-XS)/L is the panel width.  L must be at
    C        least 5 .
    C
    C     LBDCND
    C        Indicates the type of boundary conditions at X = XS and X = XF.
    C
    C        = 0  If the solution is periodic in X, i.e.
    C             U(L+I,J,K) = U(I,J,K).
    C        = 1  If the solution is specified at X = XS and X = XF.
    C        = 2  If the solution is specified at X = XS and the derivative
    C             of the solution with respect to X is specified at X = XF.
    C        = 3  If the derivative of the solution with respect to X is
    C             specified at X = XS and X = XF.
    C        = 4  If the derivative of the solution with respect to X is
    C             specified at X = XS and the solution is specified at X=XF.
    C
    C     BDXS
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to X at X = XS.
    C        when LBDCND = 3 or 4,
    C
    C             BDXS(J,K) = (d/dX)U(XS,Y(J),Z(K)), J=1,2,...,M+1,
    C                                                K=1,2,...,N+1.
    C
    C        When LBDCND has any other value, BDXS is a dummy variable.
    C        BDXS must be dimensioned at least (M+1)*(N+1).
    C
    C     BDXF
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to X at X = XF.
    C        When LBDCND = 2 or 3,
    C
    C             BDXF(J,K) = (d/dX)U(XF,Y(J),Z(K)), J=1,2,...,M+1,
    C                                                K=1,2,...,N+1.
    C
    C        When LBDCND has any other value, BDXF is a dummy variable.
    C        BDXF must be dimensioned at least (M+1)*(N+1).
    C
    C     YS,YF
    C        The range of Y, i.e. YS .LE. Y .LE. YF.
    C        YS must be less than YF.
    C
    C     M
    C        The number of panels into which the interval (YS,YF) is
    C        subdivided.  Hence, there will be M+1 grid points in the
    C        Y-direction given by Y(J) = YS+(J-1)DY for J=1,2,...,M+1,
    C        where DY = (YF-YS)/M is the panel width.  M must be at
    C        least 5 .
    C
    C     MBDCND
    C        Indicates the type of boundary conditions at Y = YS and Y = YF.
    C
    C        = 0  If the solution is periodic in Y, i.e.
    C             U(I,M+J,K) = U(I,J,K).
    C        = 1  If the solution is specified at Y = YS and Y = YF.
    C        = 2  If the solution is specified at Y = YS and the derivative
    C             of the solution with respect to Y is specified at Y = YF.
    C        = 3  If the derivative of the solution with respect to Y is
    C             specified at Y = YS and Y = YF.
    C        = 4  If the derivative of the solution with respect to Y is
    C             specified at Y = YS and the solution is specified at Y=YF.
    C
    C     BDYS
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to Y at Y = YS.
    C        When MBDCND = 3 or 4,
    C
    C             BDYS(I,K) = (d/dY)U(X(I),YS,Z(K)), I=1,2,...,L+1,
    C                                                K=1,2,...,N+1.
    C
    C        When MBDCND has any other value, BDYS is a dummy variable.
    C        BDYS must be dimensioned at least (L+1)*(N+1).
    C
    C     BDYF
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to Y at Y = YF.
    C        When MBDCND = 2 or 3,
    C
    C             BDYF(I,K) = (d/dY)U(X(I),YF,Z(K)), I=1,2,...,L+1,
    C                                                K=1,2,...,N+1.
    C
    C        When MBDCND has any other value, BDYF is a dummy variable.
    C        BDYF must be dimensioned at least (L+1)*(N+1).
    C
    C     ZS,ZF
    C        The range of Z, i.e. ZS .LE. Z .LE. ZF.
    C        ZS must be less than ZF.
    C
    C     N
    C        The number of panels into which the interval (ZS,ZF) is
    C        subdivided.  Hence, there will be N+1 grid points in the
    C        Z-direction given by Z(K) = ZS+(K-1)DZ for K=1,2,...,N+1,
    C        where DZ = (ZF-ZS)/N is the panel width.  N must be at least 5.
    C
    C     NBDCND
    C        Indicates the type of boundary conditions at Z = ZS and Z = ZF.
    C
    C        = 0  If the solution is periodic in Z, i.e.
    C             U(I,J,N+K) = U(I,J,K).
    C        = 1  If the solution is specified at Z = ZS and Z = ZF.
    C        = 2  If the solution is specified at Z = ZS and the derivative
    C             of the solution with respect to Z is specified at Z = ZF.
    C        = 3  If the derivative of the solution with respect to Z is
    C             specified at Z = ZS and Z = ZF.
    C        = 4  If the derivative of the solution with respect to Z is
    C             specified at Z = ZS and the solution is specified at Z=ZF.
    C
    C     BDZS
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to Z at Z = ZS.
    C        When NBDCND = 3 or 4,
    C
    C             BDZS(I,J) = (d/dZ)U(X(I),Y(J),ZS), I=1,2,...,L+1,
    C                                                J=1,2,...,M+1.
    C
    C        When NBDCND has any other value, BDZS is a dummy variable.
    C        BDZS must be dimensioned at least (L+1)*(M+1).
    C
    C     BDZF
    C        A two-dimensional array that specifies the values of the
    C        derivative of the solution with respect to Z at Z = ZF.
    C        When NBDCND = 2 or 3,
    C
    C             BDZF(I,J) = (d/dZ)U(X(I),Y(J),ZF), I=1,2,...,L+1,
    C                                                J=1,2,...,M+1.
    C
    C        When NBDCND has any other value, BDZF is a dummy variable.
    C        BDZF must be dimensioned at least (L+1)*(M+1).
    C
    C     ELMBDA
    C        The constant LAMBDA in the Helmholtz equation. If
    C        LAMBDA .GT. 0, a solution may not exist.  However, HW3CRT will
    C        attempt to find a solution.
    C
    C     F
    C        A three-dimensional array that specifies the values of the
    C        right side of the Helmholtz equation and boundary values (if
    C        any).  For I=2,3,...,L, J=2,3,...,M, and K=2,3,...,N
    C
    C                   F(I,J,K) = F(X(I),Y(J),Z(K)).
    C
    C        On the boundaries F is defined by
    C
    C        LBDCND      F(1,J,K)         F(L+1,J,K)
    C        ------   ---------------   ---------------
    C
    C          0      F(XS,Y(J),Z(K))   F(XS,Y(J),Z(K))
    C          1      U(XS,Y(J),Z(K))   U(XF,Y(J),Z(K))
    C          2      U(XS,Y(J),Z(K))   F(XF,Y(J),Z(K))   J=1,2,...,M+1
    C          3      F(XS,Y(J),Z(K))   F(XF,Y(J),Z(K))   K=1,2,...,N+1
    C          4      F(XS,Y(J),Z(K))   U(XF,Y(J),Z(K))
    C
    C        MBDCND      F(I,1,K)         F(I,M+1,K)
    C        ------   ---------------   ---------------
    C
    C          0      F(X(I),YS,Z(K))   F(X(I),YS,Z(K))
    C          1      U(X(I),YS,Z(K))   U(X(I),YF,Z(K))
    C          2      U(X(I),YS,Z(K))   F(X(I),YF,Z(K))   I=1,2,...,L+1
    C          3      F(X(I),YS,Z(K))   F(X(I),YF,Z(K))   K=1,2,...,N+1
    C          4      F(X(I),YS,Z(K))   U(X(I),YF,Z(K))
    C
    C        NBDCND      F(I,J,1)         F(I,J,N+1)
    C        ------   ---------------   ---------------
    C
    C          0      F(X(I),Y(J),ZS)   F(X(I),Y(J),ZS)
    C          1      U(X(I),Y(J),ZS)   U(X(I),Y(J),ZF)
    C          2      U(X(I),Y(J),ZS)   F(X(I),Y(J),ZF)   I=1,2,...,L+1
    C          3      F(X(I),Y(J),ZS)   F(X(I),Y(J),ZF)   J=1,2,...,M+1
    C          4      F(X(I),Y(J),ZS)   U(X(I),Y(J),ZF)
    C
    C        F must be dimensioned at least (L+1)*(M+1)*(N+1).
    C
    C        NOTE:
    C
    C        If the table calls for both the solution U and the right side F
    C        on a boundary, then the solution must be specified.
    C
    C     LDIMF
    C        The row (or first) dimension of the arrays F,BDYS,BDYF,BDZS,
    C        and BDZF as it appears in the program calling HW3CRT. this
    C        parameter is used to specify the variable dimension of these
    C        arrays.  LDIMF must be at least L+1.
    C
    C     MDIMF
    C        The column (or second) dimension of the array F and the row (or
    C        first) dimension of the arrays BDXS and BDXF as it appears in
    C        the program calling HW3CRT.  This parameter is used to specify
    C        the variable dimension of these arrays.
    C        MDIMF 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.  The length of W must be at least 30 + L + M + 5*N
    C        + MAX(L,M,N) + 7*(INT((L+1)/2) + INT((M+1)/2))
    C
    C
    C            * * * * * *   On Output   * * * * * *
    C
    C     F
    C        Contains the solution U(I,J,K) of the finite difference
    C        approximation for the grid point (X(I),Y(J),Z(K)) for
    C        I=1,2,...,L+1, J=1,2,...,M+1, and K=1,2,...,N+1.
    C
    C     PERTRB
    C        If a combination of periodic or derivative boundary conditions
    C        is specified for a Poisson equation (LAMBDA = 0), a solution
    C        may not exist.  PERTRB is a constant, calculated and subtracted
    C        from F, which ensures that a solution exists.  PWSCRT then
    C        computes this solution, which is a least squares solution to
    C        the original approximation.  This solution is not unique and is
    C        unnormalized.  The value of PERTRB should be small compared to
    C        the right side F.  Otherwise, a solution is obtained to an
    C        essentially different problem.  This comparison should always
    C        be made to insure that a meaningful solution has been obtained.
    C
    C     IERROR
    C        An error flag that indicates invalid input parameters.  Except
    C        for numbers 0 and 12, a solution is not attempted.
    C
    C        =  0  No error
    C        =  1  XS .GE. XF
    C        =  2  L .LT. 5
    C        =  3  LBDCND .LT. 0 .OR. LBDCND .GT. 4
    C        =  4  YS .GE. YF
    C        =  5  M .LT. 5
    C        =  6  MBDCND .LT. 0 .OR. MBDCND .GT. 4
    C        =  7  ZS .GE. ZF
    C        =  8  N .LT. 5
    C        =  9  NBDCND .LT. 0 .OR. NBDCND .GT. 4
    C        = 10  LDIMF .LT. L+1
    C        = 11  MDIMF .LT. M+1
    C        = 12  LAMBDA .GT. 0
    C
    C        Since this is the only means of indicating a possibly incorrect
    C        call to HW3CRT, the user should test IERROR after the call.
    C
    C *Long Description:
    C
    C    * * * * * * *   Program Specifications    * * * * * * * * * * * *
    C
    C     Dimension of   BDXS(MDIMF,N+1),BDXF(MDIMF,N+1),BDYS(LDIMF,N+1),
    C     Arguments      BDYF(LDIMF,N+1),BDZS(LDIMF,M+1),BDZF(LDIMF,M+1),
    C                    F(LDIMF,MDIMF,N+1),W(see argument list)
    C
    C     Latest         December 1, 1978
    C     Revision
    C
    C     Subprograms    HW3CRT,POIS3D,POS3D1,TRIDQ,RFFTI,RFFTF,RFFTF1,
    C     Required       RFFTB,RFFTB1,COSTI,COST,SINTI,SINT,COSQI,COSQF,
    C                    COSQF1,COSQB,COSQB1,SINQI,SINQF,SINQB,CFFTI,
    C                    CFFTI1,CFFTB,CFFTB1,PASSB2,PASSB3,PASSB4,PASSB,
    C                    CFFTF,CFFTF1,PASSF1,PASSF2,PASSF3,PASSF4,PASSF,
    C                    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 July 1977
    C
    C     Algorithm      This subroutine defines the finite difference
    C                    equations, incorporates boundary data, and
    C                    adjusts the right side of singular systems and
    C                    then calls POIS3D to solve the system.
    C
    C     Space          7862(decimal) = 17300(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 HW3CRT is roughly proportional
    C                    to L*M*N*(log2(L)+log2(M)+5), but also depends on
    C                    input parameters LBDCND and MBDCND.  Some typical
    C                    values are listed in the table below.
    C                       The solution process employed results in a loss
    C                    of no more than three significant digits for L,M
    C                    and N as large as 32.  More detailed information
    C                    about accuracy can be found in the documentation
    C                    for subroutine POIS3D which is the routine that
    C                    actually solves the finite difference equations.
    C
    C
    C                       L(=M=N)     LBDCND(=MBDCND=NBDCND)      T(MSECS)
    C                       -------     ----------------------      --------
    C
    C                         16                  0                    300
    C                         16                  1                    302
    C                         16                  3                    348
    C                         32                  0                   1925
    C                         32                  1                   1929
    C                         32                  3                   2109
    C
    C     Portability    American National Standards Institute FORTRAN.
    C                    The machine dependent constant PI is defined in
    C                    function PIMACH.
    C
    C     Required       COS,SIN,ATAN
    C     Resident
    C     Routines
    C
    C     Reference      NONE
    C
    C    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
    C
    C***REFERENCES  (NONE)
    C***ROUTINES CALLED  POIS3D
    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***END PROLOGUE  HW3CRT