relibc_openlibm/slatec/hw3crt.f

*DECK HW3CRT
      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
C
C
      DIMENSION       BDXS(MDIMF,*)          ,BDXF(MDIMF,*)          ,
     1                BDYS(LDIMF,*)          ,BDYF(LDIMF,*)          ,
     2                BDZS(LDIMF,*)          ,BDZF(LDIMF,*)          ,
     3                F(LDIMF,MDIMF,*)       ,W(*)
C***FIRST EXECUTABLE STATEMENT  HW3CRT
      IERROR = 0
      IF (XF .LE. XS) IERROR = 1
      IF (L .LT. 5) IERROR = 2
      IF (LBDCND.LT.0 .OR. LBDCND.GT.4) IERROR = 3
      IF (YF .LE. YS) IERROR = 4
      IF (M .LT. 5) IERROR = 5
      IF (MBDCND.LT.0 .OR. MBDCND.GT.4) IERROR = 6
      IF (ZF .LE. ZS) IERROR = 7
      IF (N .LT. 5) IERROR = 8
      IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 9
      IF (LDIMF .LT. L+1) IERROR = 10
      IF (MDIMF .LT. M+1) IERROR = 11
      IF (IERROR .NE. 0) GO TO 188
      DY = (YF-YS)/M
      TWBYDY = 2./DY
      C2 = 1./(DY**2)
      MSTART = 1
      MSTOP = M
      MP1 = M+1
      MP = MBDCND+1
      GO TO (104,101,101,102,102),MP
  101 MSTART = 2
  102 GO TO (104,104,103,103,104),MP
  103 MSTOP = MP1
  104 MUNK = MSTOP-MSTART+1
      DZ = (ZF-ZS)/N
      TWBYDZ = 2./DZ
      NP = NBDCND+1
      C3 = 1./(DZ**2)
      NP1 = N+1
      NSTART = 1
      NSTOP = N
      GO TO (108,105,105,106,106),NP
  105 NSTART = 2
  106 GO TO (108,108,107,107,108),NP
  107 NSTOP = NP1
  108 NUNK = NSTOP-NSTART+1
      LP1 = L+1
      DX = (XF-XS)/L
      C1 = 1./(DX**2)
      TWBYDX = 2./DX
      LP = LBDCND+1
      LSTART = 1
      LSTOP = L
C
C     ENTER BOUNDARY DATA FOR X-BOUNDARIES.
C
      GO TO (122,109,109,112,112),LP
  109 LSTART = 2
      DO 111 J=MSTART,MSTOP
         DO 110 K=NSTART,NSTOP
            F(2,J,K) = F(2,J,K)-C1*F(1,J,K)
  110    CONTINUE
  111 CONTINUE
      GO TO 115
  112 DO 114 J=MSTART,MSTOP
         DO 113 K=NSTART,NSTOP
            F(1,J,K) = F(1,J,K)+TWBYDX*BDXS(J,K)
  113    CONTINUE
  114 CONTINUE
  115 GO TO (122,116,119,119,116),LP
  116 DO 118 J=MSTART,MSTOP
         DO 117 K=NSTART,NSTOP
            F(L,J,K) = F(L,J,K)-C1*F(LP1,J,K)
  117    CONTINUE
  118 CONTINUE
      GO TO 122
  119 LSTOP = LP1
      DO 121 J=MSTART,MSTOP
         DO 120 K=NSTART,NSTOP
            F(LP1,J,K) = F(LP1,J,K)-TWBYDX*BDXF(J,K)
  120    CONTINUE
  121 CONTINUE
  122 LUNK = LSTOP-LSTART+1
C
C     ENTER BOUNDARY DATA FOR Y-BOUNDARIES.
C
      GO TO (136,123,123,126,126),MP
  123 DO 125 I=LSTART,LSTOP
         DO 124 K=NSTART,NSTOP
            F(I,2,K) = F(I,2,K)-C2*F(I,1,K)
  124    CONTINUE
  125 CONTINUE
      GO TO 129
  126 DO 128 I=LSTART,LSTOP
         DO 127 K=NSTART,NSTOP
            F(I,1,K) = F(I,1,K)+TWBYDY*BDYS(I,K)
  127    CONTINUE
  128 CONTINUE
  129 GO TO (136,130,133,133,130),MP
  130 DO 132 I=LSTART,LSTOP
         DO 131 K=NSTART,NSTOP
            F(I,M,K) = F(I,M,K)-C2*F(I,MP1,K)
  131    CONTINUE
  132 CONTINUE
      GO TO 136
  133 DO 135 I=LSTART,LSTOP
         DO 134 K=NSTART,NSTOP
            F(I,MP1,K) = F(I,MP1,K)-TWBYDY*BDYF(I,K)
  134    CONTINUE
  135 CONTINUE
  136 CONTINUE
C
C     ENTER BOUNDARY DATA FOR Z-BOUNDARIES.
C
      GO TO (150,137,137,140,140),NP
  137 DO 139 I=LSTART,LSTOP
         DO 138 J=MSTART,MSTOP
            F(I,J,2) = F(I,J,2)-C3*F(I,J,1)
  138    CONTINUE
  139 CONTINUE
      GO TO 143
  140 DO 142 I=LSTART,LSTOP
         DO 141 J=MSTART,MSTOP
            F(I,J,1) = F(I,J,1)+TWBYDZ*BDZS(I,J)
  141    CONTINUE
  142 CONTINUE
  143 GO TO (150,144,147,147,144),NP
  144 DO 146 I=LSTART,LSTOP
         DO 145 J=MSTART,MSTOP
            F(I,J,N) = F(I,J,N)-C3*F(I,J,NP1)
  145    CONTINUE
  146 CONTINUE
      GO TO 150
  147 DO 149 I=LSTART,LSTOP
         DO 148 J=MSTART,MSTOP
            F(I,J,NP1) = F(I,J,NP1)-TWBYDZ*BDZF(I,J)
  148    CONTINUE
  149 CONTINUE
C
C     DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
C
  150 CONTINUE
      IWB = NUNK+1
      IWC = IWB+NUNK
      IWW = IWC+NUNK
      DO 151 K=1,NUNK
         I = IWC+K-1
         W(K) = C3
         W(I) = C3
         I = IWB+K-1
         W(I) = -2.*C3+ELMBDA
  151 CONTINUE
      GO TO (155,155,153,152,152),NP
  152 W(IWC) = 2.*C3
  153 GO TO (155,155,154,154,155),NP
  154 W(IWB-1) = 2.*C3
  155 CONTINUE
      PERTRB = 0.
C
C     FOR SINGULAR PROBLEMS ADJUST DATA TO INSURE A SOLUTION WILL EXIST.
C
      GO TO (156,172,172,156,172),LP
  156 GO TO (157,172,172,157,172),MP
  157 GO TO (158,172,172,158,172),NP
  158 IF (ELMBDA) 172,160,159
  159 IERROR = 12
      GO TO 172
  160 CONTINUE
      MSTPM1 = MSTOP-1
      LSTPM1 = LSTOP-1
      NSTPM1 = NSTOP-1
      XLP = (2+LP)/3
      YLP = (2+MP)/3
      ZLP = (2+NP)/3
      S1 = 0.
      DO 164 K=2,NSTPM1
         DO 162 J=2,MSTPM1
            DO 161 I=2,LSTPM1
               S1 = S1+F(I,J,K)
  161       CONTINUE
            S1 = S1+(F(1,J,K)+F(LSTOP,J,K))/XLP
  162    CONTINUE
         S2 = 0.
         DO 163 I=2,LSTPM1
            S2 = S2+F(I,1,K)+F(I,MSTOP,K)
  163    CONTINUE
         S2 = (S2+(F(1,1,K)+F(1,MSTOP,K)+F(LSTOP,1,K)+F(LSTOP,MSTOP,K))/
     1                                                          XLP)/YLP
         S1 = S1+S2
  164 CONTINUE
      S = (F(1,1,1)+F(LSTOP,1,1)+F(1,1,NSTOP)+F(LSTOP,1,NSTOP)+
     1    F(1,MSTOP,1)+F(LSTOP,MSTOP,1)+F(1,MSTOP,NSTOP)+
     2                                   F(LSTOP,MSTOP,NSTOP))/(XLP*YLP)
      DO 166 J=2,MSTPM1
         DO 165 I=2,LSTPM1
            S = S+F(I,J,1)+F(I,J,NSTOP)
  165    CONTINUE
  166 CONTINUE
      S2 = 0.
      DO 167 I=2,LSTPM1
         S2 = S2+F(I,1,1)+F(I,1,NSTOP)+F(I,MSTOP,1)+F(I,MSTOP,NSTOP)
  167 CONTINUE
      S = S2/YLP+S
      S2 = 0.
      DO 168 J=2,MSTPM1
         S2 = S2+F(1,J,1)+F(1,J,NSTOP)+F(LSTOP,J,1)+F(LSTOP,J,NSTOP)
  168 CONTINUE
      S = S2/XLP+S
      PERTRB = (S/ZLP+S1)/((LUNK+1.-XLP)*(MUNK+1.-YLP)*
     1                                              (NUNK+1.-ZLP))
      DO 171 I=1,LUNK
         DO 170 J=1,MUNK
            DO 169 K=1,NUNK
               F(I,J,K) = F(I,J,K)-PERTRB
  169       CONTINUE
  170    CONTINUE
  171 CONTINUE
  172 CONTINUE
      NPEROD = 0
      IF (NBDCND .EQ. 0) GO TO 173
      NPEROD = 1
      W(1) = 0.
      W(IWW-1) = 0.
  173 CONTINUE
      CALL POIS3D (LBDCND,LUNK,C1,MBDCND,MUNK,C2,NPEROD,NUNK,W,W(IWB),
     1             W(IWC),LDIMF,MDIMF,F(LSTART,MSTART,NSTART),IR,W(IWW))
C
C     FILL IN SIDES FOR PERIODIC BOUNDARY CONDITIONS.
C
      IF (LP .NE. 1) GO TO 180
      IF (MP .NE. 1) GO TO 175
      DO 174 K=NSTART,NSTOP
         F(1,MP1,K) = F(1,1,K)
  174 CONTINUE
      MSTOP = MP1
  175 IF (NP .NE. 1) GO TO 177
      DO 176 J=MSTART,MSTOP
         F(1,J,NP1) = F(1,J,1)
  176 CONTINUE
      NSTOP = NP1
  177 DO 179 J=MSTART,MSTOP
         DO 178 K=NSTART,NSTOP
            F(LP1,J,K) = F(1,J,K)
  178    CONTINUE
  179 CONTINUE
  180 CONTINUE
      IF (MP .NE. 1) GO TO 185
      IF (NP .NE. 1) GO TO 182
      DO 181 I=LSTART,LSTOP
         F(I,1,NP1) = F(I,1,1)
  181 CONTINUE
      NSTOP = NP1
  182 DO 184 I=LSTART,LSTOP
         DO 183 K=NSTART,NSTOP
            F(I,MP1,K) = F(I,1,K)
  183    CONTINUE
  184 CONTINUE
  185 CONTINUE
      IF (NP .NE. 1) GO TO 188
      DO 187 I=LSTART,LSTOP
         DO 186 J=MSTART,MSTOP
            F(I,J,NP1) = F(I,J,1)
  186    CONTINUE
  187 CONTINUE
  188 CONTINUE
      RETURN
      END