cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  This is the 3D fast multipole method (FMM) code for the Yukawa 
c  potential (also called modified Helmholtz, screened Coulomb, 
c  linearized Poisson-Boltzman, etc.). It is based on the new version 
c  of FMM first introduced by Greengard and Rokhlin for the Laplace
c  equation in 1997 (see Ref. 2 below). It computes the screened 
c  Coulombic interaction between n particles.
c
c  Copyright (C) Jingfang Huang 
c
c  This program is free software; you can redistribute it and/or
c  modify it under the terms of the GNU General Public License
c  as published by the Free Software Foundation; either version 2
c  of the License, or (at your option) any later version.
c
c  This program is distributed in the hope that it will be useful,
c  but WITHOUT ANY WARRANTY; without even the implied warranty of
c  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c  GNU General Public License for more details.
c
c  License details are available in license.txt, you may also write to 
c
c    The Free Software Foundation, Inc., 
c    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
c
c  Reference:
c    1. Greengard, Leslie F.; Huang, Jingfang: A new version of the 
c       fast multipole method for screened Coulomb interactions in 
c       three dimensions. J. Comput. Phys. 180 (2002), no.2, 642--658. 
c 
c    2. Greengard, Leslie; Rokhlin, Vladimir: A new version of the 
c       fast multipole method for the Laplace equation in three dimensions.  
c       Acta numerica, 1997,  229--269, Acta Numer., 6, Cambridge Univ. 
c       Press, Cambridge, 1997. 
c
c  For suggestions, comments, and bug reports, please contact 
c
c    Jingfang Huang
c    CB# 3250, Phillips Hall
c    Department of Mathematics, UNC
c    Chapel Hill, NC 27599-3250, USA.
c    Email: huang@amath.unc.edu
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE ASSIGN(LADDER,NLEV,ZAT,IOFFST,IBOX,IATADR,
     1                   NATOMS,ICNT,ICNT2)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    assigns particles to boxes.
c
c  on input:
c
c    zat: arrays of particle coordinates.
c    nlev: refinement level
c    ladder: bookkeeping array where ladder(i) is number of boxes in
c            quad-tree of depth i.
c    icnt,icnt2 are used as workspace arrays.
c    natoms : total number of particles.
c
c  on output:
c    iatadr: an ordered list of addresses of atoms.
c            the first group of addresses correspond to atoms which
c            lie in the first box, etc.
c    ioffst(i): indicates where in the array iatadr the
c               listing of particles contained in box i begins.
c    ibox(i): address of box containing the ith atom.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
      INTEGER *4 LADDER(1),NLEV,IOFFST(1),IBOX(1),IATADR(1),NATOMS
      INTEGER *4 ICNT(1),ICNT2(1)
      REAL *8 ZAT(3,1)
c
c-----local variables
c
      INTEGER *4 NBOXES,ISTART,J,IXH,IYH,IZH
      INTEGER *4 ICOL,IROW,IPLA,IADR,IDIM,IFACE
      REAL *8 H,X,Y,Z
c
c-----initialize counting arrays
c
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      ISTART = LADDER(NLEV)
      NBOXES = LADDER(NLEV+1) - ISTART
c
      DO 100 J = 1,NBOXES
        ICNT2(J) = 1
        ICNT(J) = 0
100   CONTINUE
c
      H = 1.0D0/IDIM
c
c-----find box in which jth particle lies and increment counter array
c
      DO 200 J = 1, NATOMS
        X = ZAT(1,J) + 0.5D0
        Y = ZAT(2,J) + 0.5D0
        Z = ZAT(3,J) + 0.5D0
        IXH = X/H
        IYH = Y/H
        IZH = Z/H
        IF (IXH .GE. IDIM) IXH = IDIM-1
        IF (IYH .GE. IDIM) IYH = IDIM-1
        IF (IZH .GE. IDIM) IZH = IDIM-1
        IF (IXH .LT. 0) IXH = 0
        IF (IYH .LT. 0) IYH = 0
        IF (IZH .LT. 0) IZH = 0
        ICOL = 1 + IXH
        IROW = 1 + IYH
        IPLA = 1 + IZH
        IADR = (IPLA-1)*IFACE + (IROW-1)*IDIM + ICOL
        ICNT(IADR) = ICNT(IADR) + 1
        IBOX(J) = IADR
200   CONTINUE
c
c-----compute the array ioffst.
c
      IOFFST(1) = 0
      DO 300 J = 2,NBOXES+1
        IOFFST(J) = IOFFST(J-1) + ICNT(J-1)
300   CONTINUE
c
c-----reorder addresses of atoms in array iatadr
c
      DO 400 J = 1,NATOMS
        IADR = IBOX(J)
        IATADR(IOFFST(IADR) + ICNT2(IADR)) = J
        ICNT2(IADR) = ICNT2(IADR)+1
400   CONTINUE
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE PSBKIN(NLEV,LADDER,CENTER)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine creates the data structures needed by the
c    nonadaptive fmm.
c
c  on input:
c
c    nlev: number of levels of refinement.
c
c  on output:
c
c    ladder(nlev+1): allows us to access the information pertinent to a
c	        specific level. ladder(i)+1 is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c
c     center(3,*): coordinates of box centers. more precisely, center(i,j)
c               holds the ith component of the center of box number j.
c
c     note: for a given refinement level, boxes are numbered
c     beginning with the lowest horizontal plane (z = const.)
c     and moving up. on each plane, numbering begins in the fourth
c     quadrant corner, and proceeds in the natural ordering, one
c     row at a time.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
      INTEGER *4 NLEV, LADDER(1)
      REAL *8 CENTER(3,1)
c
c-----local variables:
c
      INTEGER *4 ICOL, IROW, IPLA
      INTEGER *4 ISTART, IFACE, ILENGTH, IDIM, JJ, I, J
      REAL *8 ZERO, H, HH, HALF
c
      DATA ZERO/0.0D0/
      DATA HALF/0.5D0/
c
c-----set hh, the length of a box side.
c
      HH = 1.0D0
      H = HH
c
c-----create the array ladder
c
      J = 1
      LADDER(1) = 0
      LADDER(2) = 1
      DO I = 3, NLEV+1
        J = J*8
        LADDER(I) = LADDER(I-1) + J
      ENDDO
c
c-----create the arrays center, iparen.
c       first do for coarsest level
c
      CENTER(1,1) = ZERO
      CENTER(2,1) = ZERO
      CENTER(3,1) = ZERO
c
c-----loop over all finer levels
c
c     idim will denote number of boxes per side at the ith level
c     iface will denote number of boxes per face.
c
      IDIM = 1
      DO 4800 I = 2, NLEV
        H = H/2
        ISTART = LADDER(I)
        ILENGTH = LADDER(I+1) - ISTART
        IDIM = IDIM*2
        IFACE = IDIM*IDIM
        DO 4600 J = 1, ILENGTH
          IPLA  = 1 + (J-1)/IFACE
          JJ = J - (IPLA-1)*IFACE
          ICOL = 1 + MOD(JJ-1,IDIM )
          IROW = 1 + (JJ-1)/IDIM
          CENTER(1,ISTART+J) = -HALF + ICOL*H - H/2
          CENTER(2,ISTART+J) = -HALF + IROW*H - H/2
          CENTER(3,ISTART+J) = -HALF + IPLA*H - H/2
 4600   CONTINUE
 4800 CONTINUE
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKNBOR(IBOX, NBORS, NNBORS, IFLAG, NLEV, LADDER,
     1  CENTER, ZCNTRS)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    this subroutine creates the list of near neighbors for box ibox.
c
c  on input:
c
c    nlev: current level.
c    ibo: box under consideration, counting from 1 at level nlev.
c    iflag: selects free-space or periodic boundary conditions.
c	  iflag = 0 for free-space problems
c	  iflag = 1 for periodic problems.
c    ladder(*): box indexing array (see psbkin).
c    center(3,*): array of box center coordinates.
c
c  on output:
c
c    nbors(*): array of near neighbor addresses counting from 1 at
c              current level.
c    nnbors: number of neighbors.
c    zcntrs(3,*): coordinates of center of ith neighbor, needed
c		  for periodic calculations.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
      INTEGER *4 IBOX,NBORS(1),NNBORS,IFLAG,NLEV,LADDER(1)
      REAL *8 CENTER(3,1),ZCNTRS(3,1)
c
c-----local variables
c
      INTEGER *4 ISTART,IDIM,IFACE,I,J,K,IXH,IYH,IZH
      INTEGER *4 ICOL,IROW,IPLA,II,JJ,KK
      REAL *8 TESTX,TESTY,TESTZ,ZERO,HALFHH,CX,CY,CZ
      REAL *8 HH, H
c
      DATA ZERO/0.0D0/
c
c-----set hh, the length of a box side, nnbors, idim, iface
c
      NNBORS = 0
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      ISTART = LADDER(NLEV)
      CX = CENTER(1,ISTART+IBOX)
      CY = CENTER(2,ISTART+IBOX)
      CZ = CENTER(3,ISTART+IBOX)
      HH = 1.0D0
      HALFHH = HH/2
      H = HH/IDIM
c
c-----compute plane, row and col numbers of current box.
c
      IXH = (CX + 0.5D0)/H
      IYH = (CY + 0.5D0)/H
      IZH = (CZ + 0.5D0)/H
      ICOL = 1 + IXH
      IROW = 1 + IYH
      IPLA = 1 + IZH
c
c-----compute nbors list (nearest shell only).
c
      IF (IFLAG .EQ. 0) THEN
c
c-------free space case
c
        DO K= -1,1
          DO J= -1,1
            DO I= -1,1
c
c-------------test to see whether computed neighbor is within the
c               computational domain.
c
              NNBORS = NNBORS + 1
              TESTX = CX + I*H
              TESTY = CY + J*H
              TESTZ = CZ + K*H
              IF (DABS(TESTX) .GT. HALFHH) THEN
                NBORS(NNBORS) = 0
              ELSEIF (DABS(TESTY) .GT. HALFHH) THEN
                NBORS(NNBORS) = 0
              ELSEIF (DABS(TESTZ) .GT. HALFHH) THEN
                NBORS(NNBORS) = 0
              ELSE
                NBORS(NNBORS) = IBOX + K*IFACE + J*IDIM + I
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ELSE IF (IFLAG .EQ. 1) THEN
c
c-------periodic case:
c
        IF (NLEV .EQ. 1) THEN
          DO K= -1,1
            DO J= -1,1
              DO I= -1,1
                NNBORS = NNBORS + 1
                ZCNTRS(1,NNBORS) = I*HH
                ZCNTRS(2,NNBORS) = J*HH
                ZCNTRS(3,NNBORS) = K*HH
                NBORS(NNBORS) = 1
              ENDDO
            ENDDO
          ENDDO
          RETURN
        ELSE
          DO 104 K= -1,1
            DO 105 J= -1,1
              DO 106 I= -1,1
                NNBORS = NNBORS + 1
                TESTX = CX + I*H
                IF (TESTX .GT. HALFHH) THEN
                  ZCNTRS(1,NNBORS) = HH
                  II = ICOL + I - IDIM
                ELSEIF (TESTX .LT. -HALFHH) THEN
                  ZCNTRS(1,NNBORS) = -HH
                  II = ICOL + I + IDIM
                ELSE
                  II = ICOL + I
                  ZCNTRS(1,NNBORS) = 0.0D0
                ENDIF
c
                TESTY = CY + J*H
                IF (TESTY .GT. HALFHH) THEN
                  ZCNTRS(2,NNBORS) = HH
                  JJ = IROW + J - IDIM
                ELSEIF (TESTY .LT. -HALFHH) THEN
                  ZCNTRS(2,NNBORS) = -HH
                  JJ = IROW + J + IDIM
                ELSE
                  JJ = IROW + J
                  ZCNTRS(2,NNBORS) = 0.0D0
                ENDIF
c
                TESTZ = CZ + K*H
                IF (TESTZ .GT. HALFHH) THEN
                  ZCNTRS(3,NNBORS) = HH
                  KK = IPLA + K - IDIM
                ELSEIF (TESTZ .LT. -HALFHH) THEN
                  ZCNTRS(3,NNBORS) = -HH
                  KK = IPLA + K + IDIM
                ELSE
                  KK = IPLA + K
                  ZCNTRS(3,NNBORS) = 0.0D0
                ENDIF
c
                NBORS(NNBORS) = (KK-1)*IFACE + (JJ-1)*IDIM + II
106           CONTINUE
105         CONTINUE
104       CONTINUE
        ENDIF
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKCHILD(NLEV,LADDER,MYBOX,ICHILD)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    compute addresses of child boxes for box j.
c      (used in nonadaptive fmm.)
c
c  on input:
c
c    nlev: number of levels of refinement.
c    ladder(*): allows us to access the information
c               pertinent to a specific level.
c	        ladder(i)+1:  is the address of the first box in the
c	        unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c    ichild(8): ichild(k) is kth child in standard
c               ordering (including appropriate offset).
c
c  note: for a given refinement level, boxes are numbered
c     beginning with the lowest horizontal plane (z = const.)
c     and moving up. on each plane, numbering begins in the fourth
c     quadrant corner, and proceeds in the natural ordering, one
c     row at a time.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),ICHILD(8),MYBOX
c
c-----local variables:
c
      INTEGER *4 IDIM,IFACE,IPLA,JJ,IPLAC
      INTEGER *4 JSTART,JBOX,ICOL,IROW,K,ISTART
      INTEGER *4 IC,JC,KC,IROWC,ICOLC
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
c-----create ichild array in standard ordering.
c
      K = 0
      ISTART = LADDER(NLEV+1)
      DO IC = 0,1
        IPLAC = IPLA*2 - 1 + IC
        DO JC = 0,1
          IROWC = IROW*2 - 1 + JC
          DO KC = 0,1
            ICOLC = ICOL*2 - 1 + KC
            K = K+1
            ICHILD(K) = ISTART + (IPLAC-1)*IFACE*4 +
     1        (IROWC-1)*IDIM*2 + ICOLC
          ENDDO
        ENDDO
      ENDDO
c
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKUPLIST(NLEV,LADDER,MYBOX,IUALL,NUALL,IXUALL,IYUALL,
     1  IU1234,NU1234,IX1234,IY1234)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     compute iuall and iu1234 interaction lists.
c
c  on input:
c
c     nlev: refinement level of box mybox (one higher than uplists.)
c     ladder(*): allows us to access the information pertinent to a
c	         specific level. ladder(i)+1:  is the address of the
c                first box in the unit cell at the ith refinement level.
c                ladder(i+1) - ladder(i): is the number of boxes in
c                the unit cell at the ith refinement level.
c     mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c     iuall(*) (integer)
c               boxes at child level receiving up expansion from all
c               eight source box children.
c     iu1234(*) (integer)
c               boxes at child level receiving up expansion from lower
c               set of source box children (i.e. 1,2,3,4).
c     nuall     (integer)
c               number of boxes in iuall list.
c     nu1234    (integer)
c               number of boxes in iu1234 lists.
c     ixuall(*)  (integer)
c               integer x offset of target boxes in iuall.
c     iyuall(*)  (integer)
c               integer y offset of target boxes in iuall.
c     ix1234(*)  (integer)
c               integer x offset of target boxes in iu1234.
c     iy1234(*)  (integer)
c               integer y offset of target boxes in iu1234.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 IUALL(36),NUALL,IXUALL(36),IYUALL(36)
      INTEGER *4 IU1234(16),NU1234,IX1234(16),IY1234(16)
c
c-----local variables:
c
      INTEGER *4 ISTART,JBOX,ICOL,IROW
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 IPLAC,IROW1,ICOL1,JMIN,JMAX,KMIN,KMAX
      INTEGER *4 IROWC,ICOLC
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NUALL = 0
      NU1234 = 0
      IF (IPLA.EQ.IDIM) THEN
        RETURN
      ENDIF
c
c-----compute iuall list.
c
      ISTART = LADDER(NLEV+1)
      IDIM = 2**NLEV
      IFACE = IDIM*IDIM
      IPLAC = IPLA*2 - 1 + 3
      IROW1 = IROW*2 - 1
      ICOL1 = ICOL*2 - 1
      JMIN = MAX(1,IROW1-2)
      JMAX = MIN(IDIM,IROW1+3)
      KMIN = MAX(1,ICOL1-2)
      KMAX = MIN(IDIM,ICOL1+3)
      DO IROWC = JMIN,JMAX
        DO ICOLC = KMIN,KMAX
          NUALL = NUALL+1
          IUALL(NUALL) = ISTART + (IPLAC-1)*IFACE +
     1      (IROWC-1)*IDIM + ICOLC
          IXUALL(NUALL) = ICOLC - ICOL1
          IYUALL(NUALL) = IROWC - IROW1
        ENDDO
      ENDDO
c
c-----compute iu1234 list.
c
      JMIN = MAX(1,IROW1-1)
      JMAX = MIN(IDIM,IROW1+2)
      KMIN = MAX(1,ICOL1-1)
      KMAX = MIN(IDIM,ICOL1+2)
      DO IROWC = JMIN,JMAX
        DO ICOLC = KMIN,KMAX
          NU1234 = NU1234+1
          IU1234(NU1234) = ISTART + (IPLAC-2)*IFACE +
     1      (IROWC-1)*IDIM + ICOLC
          IX1234(NU1234) = ICOLC - ICOL1
          IY1234(NU1234) = IROWC - IROW1
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKDNLIST(NLEV,LADDER,MYBOX,IDALL,NDALL,IXDALL,IYDALL,
     1                ID5678,ND5678,IX5678,IY5678)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c    compute idall and id5678 interaction lists.
c      (used in nonadaptive fmm.)
c
c  on input:
c
c    nlev: refinement level of box mybox (one higher than uplist.)
c    ladder(*): allows us to access the information pertinent to a
c               specific level. ladder(i)+1:  is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c     idall(*) (integer)
c               boxes at child level receiving down expansion from all
c               eight source box children.
c     id5678(*) (integer)
c               boxes at child level receiving down expansion from
c               upper set of source box children (i.e. 5,6,7,8).
c     ndall    (integer)
c               number of boxes in idnall list.
c     nd5678    (integer)
c               number of boxes in id5678 list.
c     ixdall(*)  (integer)
c               integer x offset of target boxes in idnall.
c     iydall(*)  (integer)
c               integer y offset of target boxes in idnall.
c     ix5678(*)  (integer)
c               integer x offset of target boxes in id5678.
c     iy5678(*)  (integer)
c               integer y offset of target boxes in id5678.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 IDALL(36),NDALL,IXDALL(36),IYDALL(36)
      INTEGER *4 ID5678(16),ND5678,IX5678(16),IY5678(16)
c
c-----local variables.
c
      INTEGER *4 ISTART,JBOX,ICOL,IROW
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 IPLAC,IROW1,ICOL1,JMIN,JMAX,KMIN,KMAX
      INTEGER *4 IROWC,ICOLC
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NDALL = 0
      ND5678 = 0
      IF (IPLA.EQ.1) THEN
        RETURN
      ENDIF
c
c-----compute idall list.
c
      ISTART = LADDER(NLEV+1)
      IDIM = 2**NLEV
      IFACE = IDIM*IDIM
      IPLAC = IPLA*2 - 3
      IROW1 = IROW*2 - 1
      ICOL1 = ICOL*2 - 1
      JMIN = MAX(1,IROW1-2)
      JMAX = MIN(IDIM,IROW1+3)
      KMIN = MAX(1,ICOL1-2)
      KMAX = MIN(IDIM,ICOL1+3)
      DO IROWC = JMIN,JMAX
        DO ICOLC = KMIN,KMAX
          NDALL = NDALL+1
          IDALL(NDALL) = ISTART + (IPLAC-1)*IFACE +
     1      (IROWC-1)*IDIM + ICOLC
          IXDALL(NDALL) = ICOLC - ICOL1
          IYDALL(NDALL) = IROWC - IROW1
        ENDDO
      ENDDO
c
c-----compute id5678 list.
c
      JMIN = MAX(1,IROW1-1)
      JMAX = MIN(IDIM,IROW1+2)
      KMIN = MAX(1,ICOL1-1)
      KMAX = MIN(IDIM,ICOL1+2)
      DO IROWC = JMIN,JMAX
        DO ICOLC = KMIN,KMAX
          ND5678 = ND5678+1
          ID5678(ND5678) = ISTART + (IPLAC)*IFACE +
     1      (IROWC-1)*IDIM + ICOLC
          IX5678(ND5678) = ICOLC - ICOL1
          IY5678(ND5678) = IROWC - IROW1
        ENDDO
      ENDDO
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKNOLIST(NLEV,LADDER,MYBOX,INALL,NNALL,IXNALL,IYNALL,
     1  IN1256,NN1256,IX1256,IY1256,IN12,NN12,IX12,IY12,
     2  IN56,NN56,IX56,IY56)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     compute north interaction lists.
c
c  on input:
c
c    nlev: refinement level of box mybox (one higher than uplist.)
c    ladder(*): allows us to access the information pertinent to a
c               specific level. ladder(i)+1:  is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c     inall(*) (integer)
c               boxes at child level receiving north expansion from all
c               eight source box children.
c     nnall    (integer)
c               number of boxes in inall list.
c     ixnall(*) (integer)
c               integer x offset of target boxes in inall.
c     iynall(*) (integer)
c               integer y offset of target boxes in inall.
c
c     likewise for all other lists (see algorithm description)...
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 INALL(24),NNALL,IXNALL(24),IYNALL(24)
      INTEGER *4 IN1256(8),NN1256,IX1256(8),IY1256(8)
      INTEGER *4 IN12(4),NN12,IX12(4),IY12(4)
      INTEGER *4 IN56(4),NN56,IX56(4),IY56(4)
c
c-----local variables.
c
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 ISTART,JBOX,ICOL,IROW,IPLA1,IROWC,ICOL1
      INTEGER *4 JMIN,JMAX,KMIN,KMAX,IPLAC,ICOLC
c
c-----local variables:
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NNALL = 0
      NN1256 = 0
      NN12 = 0
      NN56 = 0
      IF (IROW.NE.IDIM) THEN
c
c-------compute inall list.
c
        ISTART = LADDER(NLEV+1)
        IDIM = 2**NLEV
        IFACE = IDIM*IDIM
        IPLA1 = IPLA*2 - 1
        IROWC = IROW*2 - 1 + 3
        ICOL1 = ICOL*2 - 1
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+2)
        KMIN = MAX(1,ICOL1-2)
        KMAX = MIN(IDIM,ICOL1+3)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NNALL = NNALL+1
            INALL(NNALL) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IXNALL(NNALL) = IPLAC - IPLA1
              IYNALL(NNALL) = ICOLC - ICOL1
          ENDDO
        ENDDO
c
c-------compute in1256 list.
c
        IROWC = IROW*2 - 1 + 2
        JMIN = MAX(1,IPLA1)
        JMAX = MIN(IDIM,IPLA1+1)
        KMIN = MAX(1,ICOL1-1)
        KMAX = MIN(IDIM,ICOL1+2)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NN1256 = NN1256+1
            IN1256(NN1256) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX1256(NN1256) = IPLAC - IPLA1
            IY1256(NN1256) = ICOLC - ICOL1
          ENDDO
        ENDDO
c
c-------compute in12 and in56 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1-1)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NN12 = NN12+1
            IN12(NN12) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX12(NN12) = IPLAC - IPLA1
            IY12(NN12) = ICOLC - ICOL1
          ENDDO
        ENDDO
        JMIN = MAX(1,IPLA1+2)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NN56 = NN56+1
            IN56(NN56) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX56(NN56) = IPLAC - IPLA1
            IY56(NN56) = ICOLC - ICOL1
          ENDDO
        ENDDO
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKSOLIST(NLEV,LADDER,MYBOX,ISALL,NSALL,IXSALL,IYSALL,
     3  IS3478,NS3478,IX3478,IY3478,IS34,NS34,IX34,IY34,
     4  IS78,NS78,IX78,IY78)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     compute south interaction lists.
c
c  on input:
c
c    nlev: refinement level of box mybox (one higher than uplist.)
c    ladder(*): allows us to access the information pertinent to a
c               specific level. ladder(i)+1:  is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c    isall(*) (integer)
c             boxes at child level receiving south expansion from all
c             eight source box children.
c    nsall    (integer)
c             number of boxes in isall list.
c    ixsall(*) (integer)
c              integer x offset of target boxes in isall.
c    iysall(*) (integer)
c              integer y offset of target boxes in isall.
c
c    likewise for all other lists (see algorithm description)...
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 ISALL(24),NSALL,IXSALL(24),IYSALL(24)
      INTEGER *4 IS3478(8),NS3478,IX3478(8),IY3478(8)
      INTEGER *4 IS34(4),NS34,IX34(4),IY34(4)
      INTEGER *4 IS78(4),NS78,IX78(4),IY78(4)
c
c-----local variables.
c
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 ISTART,JBOX,ICOL,IROW,IPLA1,IROWC,ICOL1
      INTEGER *4 JMIN,JMAX,KMIN,KMAX,IPLAC,ICOLC
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NSALL = 0
      NS3478 = 0
      NS34 = 0
      NS78 = 0
      IF (IROW.NE.1) THEN
c
c-------compute isall list.
c
        ISTART = LADDER(NLEV+1)
        IDIM = 2**NLEV
        IFACE = IDIM*IDIM
        IPLA1 = IPLA*2 - 1
        IROWC = IROW*2 - 1 - 2
        ICOL1 = ICOL*2 - 1
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+2)
        KMIN = MAX(1,ICOL1-2)
        KMAX = MIN(IDIM,ICOL1+3)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NSALL = NSALL+1
            ISALL(NSALL) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IXSALL(NSALL) = IPLAC - IPLA1
            IYSALL(NSALL) = ICOLC - ICOL1
          ENDDO
        ENDDO
c
c-------compute is3478 list.
c
        IROWC = IROW*2 - 1 - 1
        JMIN = MAX(1,IPLA1)
        JMAX = MIN(IDIM,IPLA1+1)
        KMIN = MAX(1,ICOL1-1)
        KMAX = MIN(IDIM,ICOL1+2)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NS3478 = NS3478+1
            IS3478(NS3478) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX3478(NS3478) = IPLAC - IPLA1
            IY3478(NS3478) = ICOLC - ICOL1
          ENDDO
        ENDDO
c
c-------compute is34 and is78 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1-1)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NS34 = NS34+1
            IS34(NS34) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX34(NS34) = IPLAC - IPLA1
            IY34(NS34) = ICOLC - ICOL1
          ENDDO
        ENDDO
        JMIN = MAX(1,IPLA1+2)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          DO ICOLC = KMIN,KMAX
            NS78 = NS78+1
            IS78(NS78) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX78(NS78) = IPLAC - IPLA1
            IY78(NS78) = ICOLC - ICOL1
          ENDDO
        ENDDO
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKEALIST(NLEV,LADDER,MYBOX,IEALL,NEALL,IXEALL,IYEALL,
     1  IE1357,NE1357,IX1357,IY1357,IE13,NE13,IX13,IY13,
     2  IE57,NE57,IX57,IY57,IE1,NE1,IX1,IY1,IE3,NE3,IX3,IY3,
     3  IE5,NE5,IX5,IY5,IE7,NE7,IX7,IY7)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     compute east interaction lists.
c
c  on input:
c
c    nlev: refinement level of box mybox (one higher than uplist.)
c    ladder(*): allows us to access the information pertinent to a
c               specific level. ladder(i)+1:  is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c     ieall(*) (integer)
c               boxes at child level receiving east expansion from all
c               eight source box children.
c     neall    (integer)
c               number of boxes in ieall list.
c     ixeall(*) (integer)
c               integer x offset of target boxes in ieall list.
c     iyeall(*) (integer)
c               integer y offset of target boxes in ieall list.
c
c     likewise for all other lists (see algorithm description)...
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 IEALL(16),NEALL,IXEALL(16),IYEALL(16)
      INTEGER *4 IE1357(4),NE1357,IX1357(4),IY1357(4)
      INTEGER *4 IE13(2),NE13,IX13(2),IY13(2)
      INTEGER *4 IE57(2),NE57,IX57(2),IY57(2)
      INTEGER *4 IE1(3),NE1,IX1(3),IY1(3)
      INTEGER *4 IE3(3),NE3,IX3(3),IY3(3)
      INTEGER *4 IE5(3),NE5,IX5(3),IY5(3)
      INTEGER *4 IE7(3),NE7,IX7(3),IY7(3)
c
c-----local variables.
c
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 ISTART,JBOX,ICOL,IROW,IPLA1,IROWC,ICOL1
      INTEGER *4 JMIN,JMAX,KMIN,KMAX,IPLAC,ICOLC,IROW1
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NEALL = 0
      NE1357 = 0
      NE13 = 0
      NE57 = 0
      NE1 = 0
      NE3 = 0
      NE5 = 0
      NE7 = 0
      IF (ICOL.NE.IDIM) THEN
c
c-------compute ieall list.
c
        ISTART = LADDER(NLEV+1)
        IDIM = 2**NLEV
        IFACE = IDIM*IDIM
        IPLA1 = IPLA*2 - 1
        IROW1 = IROW*2 - 1
        ICOLC = ICOL*2 - 1 + 3
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+2)
        KMIN = MAX(1,IROW1-1)
        KMAX = MIN(IDIM,IROW1+2)
        DO IPLAC = JMIN,JMAX
          DO IROWC = KMIN,KMAX
            NEALL = NEALL+1
            IEALL(NEALL) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IXEALL(NEALL) = -(IPLAC - IPLA1)
            IYEALL(NEALL) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute ie1357 list.
c
        ICOLC = ICOL*2 - 1 + 2
        DO IPLAC = IPLA1,IPLA1+1
          DO IROWC = IROW1,IROW1+1
            NE1357 = NE1357+1
            IE1357(NE1357) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX1357(NE1357) = -(IPLAC - IPLA1)
            IY1357(NE1357) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute ie13 and ie57 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1-1)
        DO IPLAC = JMIN,JMAX
          DO IROWC = IROW1,IROW1+1
            NE13 = NE13+1
            IE13(NE13) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX13(NE13) = -(IPLAC - IPLA1)
            IY13(NE13) = IROWC - IROW1
          ENDDO
        ENDDO
        JMIN = MAX(1,IPLA1+2)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          DO IROWC = IROW1,IROW1+1
            NE57 = NE57+1
            IE57(NE57) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX57(NE57) = -(IPLAC - IPLA1)
            IY57(NE57) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute ie1 and ie3 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+1)
        DO IPLAC = JMIN,JMAX
          IF (IROW1.GT.1) THEN
            IROWC = IROW1-1
            NE1 = NE1+1
            IE1(NE1) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX1(NE1) = -(IPLAC - IPLA1)
            IY1(NE1) = IROWC - IROW1
          ENDIF
          IF (IROW1+2.LE.IDIM) THEN
            IROWC = IROW1+2
            NE3 = NE3+1
            IE3(NE3) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX3(NE3) = -(IPLAC - IPLA1)
            IY3(NE3) = IROWC - IROW1
          ENDIF
        ENDDO
c
c-------compute ie5 and ie7 lists.
c
        JMIN = MAX(1,IPLA1)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          IF (IROW1.GT.1) THEN
            IROWC = IROW1-1
            NE5 = NE5+1
            IE5(NE5) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX5(NE5) = -(IPLAC - IPLA1)
            IY5(NE5) = IROWC - IROW1
          ENDIF
          IF (IROW1+2.LE.IDIM) THEN
            IROWC = IROW1+2
            NE7 = NE7+1
            IE7(NE7) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX7(NE7) = -(IPLAC - IPLA1)
            IY7(NE7) = IROWC - IROW1
          ENDIF
        ENDDO
      ENDIF
c
      RETURN
      END
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
      SUBROUTINE MKWELIST(NLEV,LADDER,MYBOX,IWALL,NWALL,IXWALL,IYWALL,
     1  IW2468,NW2468,IX2468,IY2468,IW24,NW24,IX24,IY24,
     2  IW68,NW68,IX68,IY68,IW2,NW2,IX2,IY2,
     3  IW4,NW4,IX4,IY4,IW6,NW6,IX6,IY6,IW8,NW8,IX8,IY8)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccd
c
c  purpose:
c
c     compute west interaction lists.
c
c  on input:
c
c    nlev: refinement level of box mybox (one higher than uplist.)
c    ladder(*): allows us to access the information pertinent to a
c               specific level. ladder(i)+1:  is the address of the
c               first box in the unit cell at the ith refinement level.
c               ladder(i+1) - ladder(i): is the number of boxes in
c               the unit cell at the ith refinement level.
c    mybox: box under consideration (includes offset by ladder)
c
c  on output:
c
c     iwall(*) (integer)
c               boxes at child level receiving east expansion from all
c               eight source box children.
c     nwall    (integer)
c               number of boxes in iwall list.
c     ixwall(*) (integer)
c               integer x offset of target boxes in iwall list.
c     iywall(*) (integer)
c               integer y offset of target boxes in iwall list.
c
c     likewise for all other lists (see algorithm description)...
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      IMPLICIT NONE
c
      INTEGER *4 NLEV,LADDER(1),MYBOX
      INTEGER *4 IWALL(16),NWALL,IXWALL(16),IYWALL(16)
      INTEGER *4 IW2468(4),NW2468,IX2468(4),IY2468(4)
      INTEGER *4 IW24(2),NW24,IX24(2),IY24(2)
      INTEGER *4 IW68(2),NW68,IX68(2),IY68(2)
      INTEGER *4 IW2(3),NW2,IX2(3),IY2(3)
      INTEGER *4 IW4(3),NW4,IX4(3),IY4(3)
      INTEGER *4 IW6(3),NW6,IX6(3),IY6(3)
      INTEGER *4 IW8(3),NW8,IX8(3),IY8(3)
c
c-----local variables.
c
      INTEGER *4 IDIM,IFACE,IPLA,JJ
      INTEGER *4 ISTART,JBOX,ICOL,IROW,IPLA1,IROWC,ICOL1
      INTEGER *4 JMIN,JMAX,KMIN,KMAX,IPLAC,ICOLC,IROW1
c
c-----local variables:
c
      ISTART = LADDER(NLEV)
      JBOX = MYBOX - ISTART
      IDIM = 2**(NLEV-1)
      IFACE = IDIM*IDIM
      IPLA  = 1 + (JBOX-1)/IFACE
      JJ = JBOX - (IPLA-1)*IFACE
      ICOL = 1 + MOD(JJ-1,IDIM )
      IROW = 1 + (JJ-1)/IDIM
c
      NWALL = 0
      NW2468 = 0
      NW24 = 0
      NW68 = 0
      NW2 = 0
      NW4 = 0
      NW6 = 0
      NW8 = 0
      IF (ICOL.NE.1) THEN
c
c-------compute iwall list.
c
        ISTART = LADDER(NLEV+1)
        IDIM = 2**NLEV
        IFACE = IDIM*IDIM
        IPLA1 = IPLA*2 - 1
        IROW1 = IROW*2 - 1
        ICOLC = ICOL*2 - 1 - 2
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+2)
        KMIN = MAX(1,IROW1-1)
        KMAX = MIN(IDIM,IROW1+2)
        DO IPLAC = JMIN,JMAX
          DO IROWC = KMIN,KMAX
            NWALL = NWALL+1
            IWALL(NWALL) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IXWALL(NWALL) = -(IPLAC - IPLA1)
            IYWALL(NWALL) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute iw2468 list.
c
        ICOLC = ICOL*2 - 1 - 1
        DO IPLAC = IPLA1,IPLA1+1
          DO IROWC = IROW1,IROW1+1
            NW2468 = NW2468+1
            IW2468(NW2468) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX2468(NW2468) = -(IPLAC - IPLA1)
            IY2468(NW2468) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute iw24 and iw68 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1-1)
        DO IPLAC = JMIN,JMAX
          DO IROWC = IROW1,IROW1+1
            NW24 = NW24+1
            IW24(NW24) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX24(NW24) = -(IPLAC - IPLA1)
            IY24(NW24) = IROWC - IROW1
          ENDDO
        ENDDO
        JMIN = MAX(1,IPLA1+2)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          DO IROWC = IROW1,IROW1+1
            NW68 = NW68+1
            IW68(NW68) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX68(NW68) = -(IPLAC - IPLA1)
            IY68(NW68) = IROWC - IROW1
          ENDDO
        ENDDO
c
c-------compute iw2 and iw4 lists.
c
        JMIN = MAX(1,IPLA1-1)
        JMAX = MIN(IDIM,IPLA1+1)
        DO IPLAC = JMIN,JMAX
          IF (IROW1.GT.1) THEN
            IROWC = IROW1-1
            NW2 = NW2+1
            IW2(NW2) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX2(NW2) = -(IPLAC - IPLA1)
            IY2(NW2) = IROWC - IROW1
          ENDIF
          IF (IROW1+2.LE.IDIM) THEN
            IROWC = IROW1+2
            NW4 = NW4+1
            IW4(NW4) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX4(NW4) = -(IPLAC - IPLA1)
            IY4(NW4) = IROWC - IROW1
          ENDIF
        ENDDO
c
c-------compute iw6 and iw8 lists.
c
        JMIN = MAX(1,IPLA1)
        JMAX = MIN(IDIM,IPLA1+2)
        DO IPLAC = JMIN,JMAX
          IF (IROW1.GT.1) THEN
            IROWC = IROW1-1
            NW6 = NW6+1
            IW6(NW6) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX6(NW6) = -(IPLAC - IPLA1)
            IY6(NW6) = IROWC - IROW1
          ENDIF
          IF (IROW1+2.LE.IDIM) THEN
            IROWC = IROW1+2
            NW8 = NW8+1
            IW8(NW8) = ISTART + (IPLAC-1)*IFACE +
     1        (IROWC-1)*IDIM + ICOLC
            IX8(NW8) = -(IPLAC - IPLA1)
            IY8(NW8) = IROWC - IROW1
          ENDIF
        ENDDO
      ENDIF
c
      RETURN
      END
