C
C $Id: pgridrfnc.F,v 1.3 1998/07/16 16:40:12 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     this file contains all the subroutines needed
C     for residue-wise, grid-based subsetting with periodic boundary
C     conditions. These subroutines are called in the driver subroutines
C     mixgrid and resgrid.
C
C     Written by Arjan van der Vaart, July-August 1997
C     .                               January 1998 (pbc)
C
C     note that in all the subroutines i refers to
C     the core (i=1), inner buffer (i=2) or outer buffer (i=3).
C     k is a counter for igsub, igsub is filled up with
C     atoms within the allowed range. kx is a counter for igsubx,
C     igsubx is filled up with atoms that are outside of the allowed
C     range (as defined by gridmin and gridmax).
C
C     The code is elaborate and long to save a lot of
C     if-statements. These if-statements verify if the residues
C     are within certain boundaries (set by gridmin and gridmax).
C
C     with i the number of cells in the x-direction,
C     .    j the number of cells in the y-direction,
C     .    k the number of cells in the z-direction for a subsystem,
C
C     a crude-force method would need 6*i*j*k if-statements, while this
C     method only needs 2(i*j + i*k + j*k) if-statements; for a cube
C     a crude-force method would need 6*n*n*n if-statements, this
C     method 6*n*n .
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC



CC--------------------------------------------------------------------------CC

      subroutine pgetrcorners(i, k, igsub, kx, igsubx)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     tedious and elaborate code to save half the
C     number of if-statements

C     i=1 -> core
C     i=2 -> inner buffer
C     i=3 -> outer buffer

      dimension igsub(*), igsubx(*)

C     locals:
      dimension gmin(3), gmax(3)

CC---------------------------------------------------------CC

      if (ibg(3,i).lt.0) then
         ibg3 = ibg(3,i)+ncell(3)
         gmin(3) = dbox(3)+gridmin(3,i)
      else
         ibg3 = ibg(3,i)
         gmin(3) = gridmin(3,i)
      endif
      if (ibg(2,i).lt.0) then
         ibg2 = ibg(2,i)+ncell(2)
         gmin(2) = dbox(2)+gridmin(2,i)
      else
         ibg2 = ibg(2,i)
         gmin(2) = gridmin(2,i)
      endif
      if (ibg(1,i).lt.0) then
         ibg1 = ibg(1,i)+ncell(1)
         gmin(1) = dbox(1)+gridmin(1,i)
      else
         ibg1 = ibg(1,i)
         gmin(1) = gridmin(1,i)
      endif

      if (ied(3,i).ge.ncell(3)) then
         ied3 = ied(3,i)-ncell(3)
         gmax(3) = gridmax(3,i)-dbox(3)
      else
         ied3 = ied(3,i)
         gmax(3) = gridmax(3,i)
      endif
      if (ied(2,i).ge.ncell(2)) then
         ied2 = ied(2,i)-ncell(2)
         gmax(2) = gridmax(2,i)-dbox(2)
      else
         ied2 = ied(2,i)
         gmax(2) = gridmax(2,i)
      endif
      if (ied(1,i).ge.ncell(1)) then
         ied1 = ied(1,i)-ncell(1)
         gmax(1) = gridmax(1,i)-dbox(1)
      else
         ied1 = ied(1,i)
         gmax(1) = gridmax(1,i)
      endif

      icell0 = ibg3*ncell2 + ibg2*ncell(1) + 1
      icell = icell0 + ibg1
      irese = igridradr(icell+1)-1
      do 40 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).lt.gmin(1)) goto 20
         if (gc(2,kres).lt.gmin(2)) goto 20
         if (gc(3,kres).lt.gmin(3)) goto 20
         do 10 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 10      continue
         goto 40
 20      do 30 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 30      continue
 40   continue

      icell = icell0 + ied1
      irese = igridradr(icell+1)-1
      do 80 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).gt.gmax(1)) goto 60
         if (gc(2,kres).lt.gmin(2)) goto 60
         if (gc(3,kres).lt.gmin(3)) goto 60
         do 50 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 50      continue
         goto 80
 60      do 70 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 70      continue
 80   continue

      icell0 = ibg3*ncell2 + ied2*ncell(1) + 1
      icell = icell0 + ibg1
      irese = igridradr(icell+1)-1
      do 120 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).lt.gmin(1)) goto 100
         if (gc(2,kres).gt.gmax(2)) goto 100
         if (gc(3,kres).lt.gmin(3)) goto 100
         do 90 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 90      continue
         goto 120
 100     do 110 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 110     continue
 120  continue

      icell = icell0 + ied1
      irese = igridradr(icell+1)-1
      do 160 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).gt.gmax(1)) goto 140
         if (gc(2,kres).gt.gmax(2)) goto 140
         if (gc(3,kres).lt.gmin(3)) goto 140
         do 130 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 130     continue
         goto 160
 140     do 150 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 150     continue
 160  continue

      icell0 = ied3*ncell2 + ibg2*ncell(1) + 1
      icell = icell0 + ibg1
      irese = igridradr(icell+1)-1
      do 200 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).lt.gmin(1)) goto 180
         if (gc(2,kres).lt.gmin(2)) goto 180
         if (gc(3,kres).gt.gmax(3)) goto 180
         do 170 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 170     continue
         goto 200
 180     do 190 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 190     continue
 200  continue

      icell = icell0 + ied1
      irese = igridradr(icell+1)-1
      do 240 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).gt.gmax(1)) goto 220
         if (gc(2,kres).lt.gmin(2)) goto 220
         if (gc(3,kres).gt.gmax(3)) goto 220
         do 210 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 210     continue
         goto 240
 220     do 230 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 230     continue
 240  continue

      icell0 = ied3*ncell2 + ied2*ncell(1) + 1
      icell = icell0 + ibg1
      irese = igridradr(icell+1)-1
      do 280 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).lt.gmin(1)) goto 260
         if (gc(2,kres).gt.gmax(2)) goto 260
         if (gc(3,kres).gt.gmax(3)) goto 260
         do 250 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 250     continue
         goto 280
 260     do 270 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 270     continue
 280  continue

      icell = icell0 + ied1
      irese = igridradr(icell+1)-1
      do 320 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         if (gc(1,kres).gt.gmax(1)) goto 300
         if (gc(2,kres).gt.gmax(2)) goto 300
         if (gc(3,kres).gt.gmax(3)) goto 300
         do 290 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 290     continue
         goto 320
 300     do 310 iat=irpnt(kres),iate
            kx = kx + 1
            igsubx(kx) = iat
 310     continue
 320  continue

      end

CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

      subroutine pgetrribs(i, k, igsub, kx, igsubx)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension igsub(*), igsubx(*)

C     tedious and elaborate code to save a lot
C     of if-statements

C     locals:
      dimension gmin(3), gmax(3)

CC--------------------------------------------------CC

      ixb = ibg(1,i)+1
      ixe = ied(1,i)-1
      iyb = ibg(2,i)+1
      iye = ied(2,i)-1
      izb = ibg(3,i)+1
      ize = ied(3,i)-1

      if (ibg(3,i).lt.0) then
         ibg3 = ibg(3,i)+ncell(3)
         gmin(3) = dbox(3)+gridmin(3,i)
      else
         ibg3 = ibg(3,i)
         gmin(3) = gridmin(3,i)
      endif
      if (ibg(2,i).lt.0) then
         ibg2 = ibg(2,i)+ncell(2)
         gmin(2) = dbox(2)+gridmin(2,i)
      else
         ibg2 = ibg(2,i)
         gmin(2) = gridmin(2,i)
      endif
      if (ibg(1,i).lt.0) then
         ibg1 = ibg(1,i)+ncell(1)
         gmin(1) = dbox(1)+gridmin(1,i)
      else
         ibg1 = ibg(1,i)
         gmin(1) = gridmin(1,i)
      endif

      if (ied(3,i).ge.ncell(3)) then
         ied3 = ied(3,i)-ncell(3)
         gmax(3) = gridmax(3,i)-dbox(3)
      else
         ied3 = ied(3,i)
         gmax(3) = gridmax(3,i)
      endif
      if (ied(2,i).ge.ncell(2)) then
         ied2 = ied(2,i)-ncell(2)
         gmax(2) = gridmax(2,i)-dbox(2)
      else
         ied2 = ied(2,i)
         gmax(2) = gridmax(2,i)
      endif
      if (ied(1,i).ge.ncell(1)) then
         ied1 = ied(1,i)-ncell(1)
         gmax(1) = gridmax(1,i)-dbox(1)
      else
         ied1 = ied(1,i)
         gmax(1) = gridmax(1,i)
      endif

C     zx-ribs
      icell0 = ibg3*ncell2 + ibg1 + 1
      do 50 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         irese = igridradr(icell+1)-1
         do 40 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).lt.gmin(1)) goto 20
            if (gc(3,kres).lt.gmin(3)) goto 20
            do 10 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 10         continue
            goto 40
 20         do 30 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
               continue
 30         continue
 40      continue
 50   continue

      icell0 = ibg3*ncell2 + ied1 + 1
      do 100 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         irese = igridradr(icell+1)-1
         do 90 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).gt.gmax(1)) goto 70
            if (gc(3,kres).lt.gmin(3)) goto 70
            do 60 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 60         continue
            goto 90
 70         do 80 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 80         continue
 90      continue
 100  continue

      icell0 = ied3*ncell2 + ibg1 + 1
      do 150 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         irese = igridradr(icell+1)-1
         do 140 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).lt.gmin(1)) goto 120
            if (gc(3,kres).gt.gmax(3)) goto 120
            do 110 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 110        continue
            goto 140
 120        do 130 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 130        continue
 140     continue
 150  continue

      icell0 = ied3*ncell2 + ied1 + 1
      do 200 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         irese = igridradr(icell+1)-1
         do 190 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).gt.gmax(1)) goto 170
            if (gc(3,kres).gt.gmax(3)) goto 170
            do 160 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 160        continue
            goto 190
 170        do 180 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 180        continue
 190     continue
 200  continue

C     zy-ribs
      icell0 = ibg3*ncell2 + ibg2*ncell(1) + 1
      do 250 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         irese = igridradr(icell+1)-1
         do 240 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(2,kres).lt.gmin(2)) goto 220
            if (gc(3,kres).lt.gmin(3)) goto 220
            do 210 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 210        continue
            goto 240
 220        do 230 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 230        continue
 240     continue
 250  continue

      icell0 = ibg3*ncell2 + ied2*ncell(1) + 1
      do 300 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         irese = igridradr(icell+1)-1
         do 290 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(2,kres).gt.gmax(2)) goto 270
            if (gc(3,kres).lt.gmin(3)) goto 270
            do 260 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 260        continue
            goto 290
 270        do 280 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 280        continue
 290     continue
 300  continue

      icell0 = ied3*ncell2 + ibg2*ncell(1) + 1
      do 350 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         irese = igridradr(icell+1)-1
         do 340 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(2,kres).lt.gmin(2)) goto 320
            if (gc(3,kres).gt.gmax(3)) goto 320
            do 310 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 310        continue
            goto 340
 320        do 330 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 330        continue
 340     continue
 350  continue

      icell0 = ied3*ncell2 + ied2*ncell(1) + 1
      do 400 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         irese = igridradr(icell+1)-1
         do 390 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(2,kres).gt.gmax(2)) goto 370
            if (gc(3,kres).gt.gmax(3)) goto 370
            do 360 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 360        continue
            goto 390
 370        do 380 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 380        continue
 390     continue
 400  continue

C     xy-ribs
      icell0 = ibg2*ncell(1) + ibg1 + 1
      do 450 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         irese = igridradr(icell+1)-1
         do 440 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).lt.gmin(1)) goto 420
            if (gc(2,kres).lt.gmin(2)) goto 420
            do 410 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 410        continue
            goto 440
 420        do 430 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 430        continue
 440     continue
 450  continue

      icell0 = ibg2*ncell(1) + ied1 + 1
      do 500 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         irese = igridradr(icell+1)-1
         do 490 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).gt.gmax(1)) goto 470
            if (gc(2,kres).lt.gmin(2)) goto 470
            do 460 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 460        continue
            goto 490
 470        do 480 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 480        continue
 490     continue
 500  continue

      icell0 = ied2*ncell(1) + ibg1 + 1
      do 550 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         irese = igridradr(icell+1)-1
         do 540 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).lt.gmin(1)) goto 520
            if (gc(2,kres).gt.gmax(2)) goto 520
            do 510 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 510        continue
            goto 540
 520        do 530 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 530        continue
 540     continue
 550  continue

      icell0 = ied2*ncell(1) + ied1 + 1
      do 600 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         irese = igridradr(icell+1)-1
         do 590 ires=igridradr(icell),irese
            kres = igridr(ires)
            iate = irpnt(kres+1)-1
            if (gc(1,kres).gt.gmax(1)) goto 570
            if (gc(2,kres).gt.gmax(2)) goto 570
            do 560 iat=irpnt(kres),iate
               k = k + 1
               igsub(k) = iat
 560        continue
            goto 590
 570        do 580 iat=irpnt(kres),iate
               kx = kx + 1
               igsubx(kx) = iat
 580        continue
 590     continue
 600  continue

      end


CC--------------------------------------------------------------CC
CC--------------------------------------------------------------CC

      subroutine pgetrplanes(i, k, igsub, kx, igsubx)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension igsub(*), igsubx(*)

C     locals:
      dimension gmin(3), gmax(3)

CC-------------------------------------------------------CC

      ibx = ibg(1,i)+1
      iex = ied(1,i)-1
      iby = ibg(2,i)+1
      iey = ied(2,i)-1
      ibz = ibg(3,i)+1
      iez = ied(3,i)-1

      if (ibg(3,i).lt.0) then
         ibg3 = ibg(3,i)+ncell(3)
         gmin(3) = dbox(3)+gridmin(3,i)
      else
         ibg3 = ibg(3,i)
         gmin(3) = gridmin(3,i)
      endif
      if (ibg(2,i).lt.0) then
         ibg2 = ibg(2,i)+ncell(2)
         gmin(2) = dbox(2)+gridmin(2,i)
      else
         ibg2 = ibg(2,i)
         gmin(2) = gridmin(2,i)
      endif
      if (ibg(1,i).lt.0) then
         ibg1 = ibg(1,i)+ncell(1)
         gmin(1) = dbox(1)+gridmin(1,i)
      else
         ibg1 = ibg(1,i)
         gmin(1) = gridmin(1,i)
      endif

      if (ied(3,i).ge.ncell(3)) then
         ied3 = ied(3,i)-ncell(3)
         gmax(3) = gridmax(3,i)-dbox(3)
      else
         ied3 = ied(3,i)
         gmax(3) = gridmax(3,i)
      endif
      if (ied(2,i).ge.ncell(2)) then
         ied2 = ied(2,i)-ncell(2)
         gmax(2) = gridmax(2,i)-dbox(2)
      else
         ied2 = ied(2,i)
         gmax(2) = gridmax(2,i)
      endif
      if (ied(1,i).ge.ncell(1)) then
         ied1 = ied(1,i)-ncell(1)
         gmax(1) = gridmax(1,i)-dbox(1)
      else
         ied1 = ied(1,i)
         gmax(1) = gridmax(1,i)
      endif

C     z planes:

      it = ibg3*ncell2 + 1
      do 60 iy=iby,iey
         iiy = mod(iy+ncell(2),ncell(2))
         icell0 = iiy*ncell(1) + it
         do 50 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            irese = igridradr(icell+1)-1
            do 40 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(3,kres).lt.gmin(3)) goto 20
               do 10 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 10            continue
               goto 40
 20            do 30 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 30            continue
 40         continue
 50      continue
 60   continue

      it = ied3*ncell2 + 1
      do 120 iy=iby,iey
         iiy = mod(iy+ncell(2),ncell(2))
         icell0 = iiy*ncell(1) + it
         do 110 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            irese = igridradr(icell+1)-1
            do 100 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(3,kres).gt.gmax(3)) goto 80
               do 70 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 70            continue
               goto 100
 80            do 90 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 90            continue
 100        continue
 110     continue
 120  continue

C     y planes:

      it = ibg2*ncell(1) + 1
      do 180 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 170 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            irese = igridradr(icell+1)-1
            do 160 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(2,kres).lt.gmin(2)) goto 140
               do 130 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 130           continue
               goto 160
 140           do 150 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 150           continue
 160        continue
 170     continue
 180  continue

      it = ied2*ncell(1) + 1
      do 240 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 230 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            irese = igridradr(icell+1)-1
            do 220 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(2,kres).gt.gmax(2)) goto 200
               do 190 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 190           continue
               goto 220
 200           do 210 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 210           continue
 220        continue
 230     continue
 240  continue

C     x planes:

      it = ibg1 + 1
      do 300 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 290 iy=iby,iey
            iiy = mod(iy+ncell(2),ncell(2))
            icell = icell0 + iiy*ncell(1)
            irese = igridradr(icell+1)-1
            do 280 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(1,kres).lt.gmin(1)) goto 260
               do 250 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 250           continue
               goto 280
 260           do 270 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 270           continue
 280        continue
 290     continue
 300  continue

      it = ied1 + 1
      do 360 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 350 iy=iby,iey
            iiy = mod(iy+ncell(2),ncell(2))
            icell = icell0 + iiy*ncell(1)
            irese = igridradr(icell+1)-1
            do 340 ires=igridradr(icell),irese
               kres = igridr(ires)
               iate = irpnt(kres+1)-1
               if (gc(1,kres).gt.gmax(1)) goto 320
               do 310 iat=irpnt(kres),iate
                  k = k + 1
                  igsub(k) = iat
 310           continue
               goto 340
 320           do 330 iat=irpnt(kres),iate
                  kx = kx + 1
                  igsubx(kx) = iat
 330           continue
 340        continue
 350     continue
 360  continue

      end


CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

      subroutine pgetrcore(k, igsub)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension igsub(*)

CC---------------------------------------------------CC

      ixb = ibg(1,1)+1
      ixe = ied(1,1)-1
      iyb = ibg(2,1)+1
      iye = ied(2,1)-1
      izb = ibg(3,1)+1
      ize = ied(3,1)-1

      do 30 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         izn = iiz*ncell2 + 1
         do 20 iy=iyb,iye
            iiy = mod(iy+ncell(2),ncell(2))
            icell0 = izn + iiy*ncell(1)
            do 10 ix=ixb,ixe
               iix = mod(ix+ncell(1),ncell(1))
               icell = icell0 + iix
               call getrat(icell, k, igsub)
 10         continue
 20      continue
 30   continue

      end

CC-----------------------------------------------------------------CC
CC-----------------------------------------------------------------CC

      subroutine pgetrat(icell, k, igsub)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension igsub(*)

      irese = igridradr(icell+1)-1
      do 20 ires=igridradr(icell),irese
         kres = igridr(ires)
         iate = irpnt(kres+1)-1
         do 10 iat=irpnt(kres),iate
            k = k + 1
            igsub(k) = iat
 10      continue
 20   continue

      end

CC----------------------------------------------------------------CC
CC----------------------------------------------------------------CC

      subroutine pgetrinner(i, k, igsub)

      implicit double precision (a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

      dimension igsub(*)


CC---------------------------------------------CC

C     z layer
      ixb = ibg(1,i)+1
      ixe = ied(1,i)-1
      iyb = ibg(2,i)+1
      iye = ied(2,i)-1
      izb1 = ibg(3,i)+1
      ize1 = ibg(3,i-1)-1
      izb2 = ied(3,i-1)+1
      ize2 = ied(3,i)-1

      do iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         iyn = iiy*ncell(1) + 1
         do ix=ixb,ixe
            iix = mod(ix+ncell(1),ncell(1))
            icell0 = iyn + iix
            do iz=izb1,ize1
               iiz = mod(iz+ncell(3),ncell(3))
               icell = icell0 + iiz*ncell2
               call getrat(icell, k, igsub)
            enddo
            do iz=izb2,ize2
               iiz = mod(iz+ncell(3),ncell(3))
               icell = icell0 + iiz*ncell2
               call getrat(icell, k, igsub)
            enddo
         enddo
      enddo

C     y layer

      izb = ibg(3,i-1)
      ize = ied(3,i-1)
      iyb1 = ibg(2,i)+1
      iye1 = ibg(2,i-1)-1
      iyb2 = ied(2,i-1)+1
      iye2 = ied(2,i)-1

      do ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         ixn = iix + 1
         do iz=izb,ize
            iiz = mod(iz+ncell(3),ncell(3))
            icell0 = ixn + iiz*ncell2
            do iy=iyb1,iye1
               iiy = mod(iy+ncell(2),ncell(2))
               icell = icell0 + iiy*ncell(1)
               call getrat(icell, k, igsub)
            enddo
            do iy=iyb2,iye2
               iiy = mod(iy+ncell(2),ncell(2))
               icell = icell0 + iiy*ncell(1)
               call getrat(icell, k, igsub)
            enddo
         enddo
      enddo

C     x layer

      iyb = ibg(2,i-1)
      iye = ied(2,i-1)
      izb = ibg(3,i-1)
      ize = ied(3,i-1)
      ixb1 = ibg(1,i)+1
      ixe1 = ibg(1,i-1)-1
      ixb2 = ied(1,i-1)+1
      ixe2 = ied(1,i)-1

      do iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         iyn = iiy*ncell(1) + 1
         do iz=izb,ize
            iiz = mod(iz+ncell(3),ncell(3))
            icell0 = iyn + iiz*ncell2
            do ix=ixb1,ixe1
               iix = mod(ix+ncell(1),ncell(1))
               icell = icell0 + iix
               call getrat(icell, k, igsub)
            enddo
            do ix=ixb2,ixe2
               iix = mod(ix+ncell(1),ncell(1))
               icell = icell0 + iix
               call getrat(icell, k, igsub)
            enddo
         enddo
      enddo

      end

CC------------------------------------------------------------------------CC