C
C $Id: pgridafnc.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 atom-wise, grid-based subsetting with periodic boundary
C     conditions. These subroutines are called in the driver subroutines
C     atgrid and mixgrid.
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 atoms
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 pgetacorners(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
      iate = igridadr(icell+1)-1
      do 20 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).lt.gmin(1)) goto 10
         if (xyz(2,kat).lt.gmin(2)) goto 10
         if (xyz(3,kat).lt.gmin(3)) goto 10
         k = k + 1
         igsub(k) = kat
         goto 20
 10      kx = kx + 1
         igsubx(kx) = kat
 20   continue

      icell = icell0 + ied1
      iate = igridadr(icell+1)-1
      do 40 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).gt.gmax(1)) goto 30
         if (xyz(2,kat).lt.gmin(2)) goto 30
         if (xyz(3,kat).lt.gmin(3)) goto 30
         k = k + 1
         igsub(k) = kat
         goto 40
 30      kx = kx + 1
         igsubx(kx) = kat
 40   continue

      icell0 = ibg3*ncell2 + ied2*ncell(1) + 1
      icell = icell0 + ibg1
      iate = igridadr(icell+1)-1
      do 60 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).lt.gmin(1)) goto 50
         if (xyz(2,kat).gt.gmax(2)) goto 50
         if (xyz(3,kat).lt.gmin(3)) goto 50
         k = k + 1
         igsub(k) = kat
         goto 60
 50      kx = kx + 1
         igsubx(kx) = kat
 60   continue

      icell = icell0 + ied1
      iate = igridadr(icell+1)-1
      do 80 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).gt.gmax(1)) goto 70
         if (xyz(2,kat).gt.gmax(2)) goto 70
         if (xyz(3,kat).lt.gmin(3)) goto 70
         k = k + 1
         igsub(k) = kat
         goto 80
 70      kx = kx + 1
         igsubx(kx) = kat
 80   continue

      icell0 = ied3*ncell2 + ibg2*ncell(1) + 1
      icell = icell0 + ibg1
      iate = igridadr(icell+1)-1
      do 100 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).lt.gmin(1)) goto 90
         if (xyz(2,kat).lt.gmin(2)) goto 90
         if (xyz(3,kat).gt.gmax(3)) goto 90
         k = k + 1
         igsub(k) = kat
         goto 100
 90      kx = kx + 1
         igsubx(kx) = kat
 100  continue

      icell = icell0 + ied1
      iate = igridadr(icell+1)-1
      do 120 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).gt.gmax(1)) goto 110
         if (xyz(2,kat).lt.gmin(2)) goto 110
         if (xyz(3,kat).gt.gmax(3)) goto 110
         k = k + 1
         igsub(k) = kat
         goto 120
 110     kx = kx + 1
         igsubx(kx) = kat
 120  continue

      icell0 = ied3*ncell2 + ied2*ncell(1) + 1
      icell = icell0 + ibg1
      iate = igridadr(icell+1)-1
      do 140 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).lt.gmin(1)) goto 130
         if (xyz(2,kat).gt.gmax(2)) goto 130
         if (xyz(3,kat).gt.gmax(3)) goto 130
         k = k + 1
         igsub(k) = kat
         goto 140
 130     kx = kx + 1
         igsubx(kx) = kat
 140  continue

      icell = icell0 + ied1
      iate = igridadr(icell+1)-1
      do 160 iat=igridadr(icell),iate
         kat = igrid(iat)
         if (xyz(1,kat).gt.gmax(1)) goto 170
         if (xyz(2,kat).gt.gmax(2)) goto 170
         if (xyz(3,kat).gt.gmax(3)) goto 170
         k = k + 1
         igsub(k) = kat
         goto 160
 170     kx = kx + 1
         igsubx(kx) = kat
 160  continue

      end

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

      subroutine pgetaribs(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 30 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         iate = igridadr(icell+1)-1
         do 20 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).lt.gmin(1)) goto 10
            if (xyz(3,kat).lt.gmin(3)) goto 10
            k = k + 1
            igsub(k) = kat
            goto 20
 10         kx = kx + 1
            igsubx(kx) = kat
 20      continue
 30   continue

      icell0 = ibg3*ncell2 + ied1 + 1
      do 60 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         iate = igridadr(icell+1)-1
         do 50 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).gt.gmax(1)) goto 40
            if (xyz(3,kat).lt.gmin(3)) goto 40
            k = k + 1
            igsub(k) = kat
            goto 50
 40         kx = kx + 1
            igsubx(kx) = kat
 50      continue
 60   continue

      icell0 = ied3*ncell2 + ibg1 + 1
      do 90 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         iate = igridadr(icell+1)-1
         do 80 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).lt.gmin(1)) goto 70
            if (xyz(3,kat).gt.gmax(3)) goto 70
            k = k + 1
            igsub(k) = kat
            goto 80
 70         kx = kx + 1
            igsubx(kx) = kat
 80      continue
 90   continue

      icell0 = ied3*ncell2 + ied1 + 1
      do 120 iy=iyb,iye
         iiy = mod(iy+ncell(2),ncell(2))
         icell = icell0 + iiy*ncell(1)
         iate = igridadr(icell+1)-1
         do 110 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).gt.gmax(1)) goto 100
            if (xyz(3,kat).gt.gmax(3)) goto 100
            k = k + 1
            igsub(k) = kat
            goto 110
 100        kx = kx + 1
            igsubx(kx) = kat
 110     continue
 120  continue

C     zy-ribs

      icell0 = ibg3*ncell2 + ibg2*ncell(1) + 1
      do 150 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         iate = igridadr(icell+1)-1
         do 140 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(2,kat).lt.gmin(2)) goto 130
            if (xyz(3,kat).lt.gmin(3)) goto 130
            k = k + 1
            igsub(k) = kat
            goto 140
 130        kx = kx + 1
            igsubx(kx) = kat
 140     continue
 150  continue

      icell0 = ibg3*ncell2 + ied2*ncell(1) + 1
      do 180 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         iate = igridadr(icell+1)-1
         do 170 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(2,kat).gt.gmax(2)) goto 160
            if (xyz(3,kat).lt.gmin(3)) goto 160
            k = k + 1
            igsub(k) = kat
            goto 170
 160        kx = kx + 1
            igsubx(kx) = kat
 170     continue
 180  continue

      icell0 = ied3*ncell2 + ibg2*ncell(1) + 1
      do 210 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         iate = igridadr(icell+1)-1
         do 200 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(2,kat).lt.gmin(2)) goto 190
            if (xyz(3,kat).gt.gmax(3)) goto 190
            k = k + 1
            igsub(k) = kat
            goto 200
 190        kx = kx + 1
            igsubx(kx) = kat
 200     continue
 210  continue

      icell0 = ied3*ncell2 + ied2*ncell(1) + 1
      do 240 ix=ixb,ixe
         iix = mod(ix+ncell(1),ncell(1))
         icell = icell0 + iix
         iate = igridadr(icell+1)-1
         do 230 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(2,kat).gt.gmax(2)) goto 220
            if (xyz(3,kat).gt.gmax(3)) goto 220
            k = k + 1
            igsub(k) = kat
            goto 230
 220        kx = kx + 1
            igsubx(kx) = kat
 230     continue
 240  continue

C     xy-ribs

      icell0 = ibg2*ncell(1) + ibg1 + 1
      do 270 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         iate = igridadr(icell+1)-1
         do 260 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).lt.gmin(1)) goto 250
            if (xyz(2,kat).lt.gmin(2)) goto 250
            k = k + 1
            igsub(k) = kat
            goto 260
 250        kx = kx + 1
            igsubx(kx) = kat
 260     continue
 270  continue

      icell0 = ibg2*ncell(1) + ied1 + 1
      do 300 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         iate = igridadr(icell+1)-1
         do 290 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).gt.gmax(1)) goto 280
            if (xyz(2,kat).lt.gmin(2)) goto 280
            k = k + 1
            igsub(k) = kat
            goto 290
 280        kx = kx + 1
            igsubx(kx) = kat
 290     continue
 300  continue

      icell0 = ied2*ncell(1) + ibg1 + 1
      do 330 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         iate = igridadr(icell+1)-1
         do 320 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).lt.gmin(1)) goto 310
            if (xyz(2,kat).gt.gmax(2)) goto 310
            k = k + 1
            igsub(k) = kat
            goto 320
 310        kx = kx + 1
            igsubx(kx) = kat
 320     continue
 330  continue

      icell0 = ied2*ncell(1) + ied1 + 1
      do 360 iz=izb,ize
         iiz = mod(iz+ncell(3),ncell(3))
         icell = icell0 + iiz*ncell2
         iate = igridadr(icell+1)-1
         do 350 iat=igridadr(icell),iate
            kat = igrid(iat)
            if (xyz(1,kat).gt.gmax(1)) goto 340
            if (xyz(2,kat).gt.gmax(2)) goto 340
            k = k + 1
            igsub(k) = kat
            goto 350
 340        kx = kx + 1
            igsubx(kx) = kat
 350     continue
 360  continue

      end

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

      subroutine pgetaplanes(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 40 iy=iby,iey
         iiy = mod(iy+ncell(2),ncell(2))
         icell0 = iiy*ncell(1) + it
         do 30 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            iate = igridadr(icell+1)-1
            do 20 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(3,kat).lt.gmin(3)) goto 10
               k = k + 1
               igsub(k) = kat
               goto 20
 10            kx = kx + 1
               igsubx(kx) = kat
 20         continue
 30      continue
 40   continue

      it = ied3*ncell2 + 1
      do 80 iy=iby,iey
         iiy = mod(iy+ncell(2),ncell(2))
         icell0 = iiy*ncell(1) + it
         do 70 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            iate = igridadr(icell+1)-1
            do 60 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(3,kat).gt.gmax(3)) goto 50
               k = k + 1
               igsub(k) = kat
               goto 60
 50            kx = kx + 1
               igsubx(kx) = kat
 60         continue
 70      continue
 80   continue

C     y planes:

      it = ibg2*ncell(1) + 1
      do 120 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 110 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            iate = igridadr(icell+1)-1
            do 100 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).lt.gmin(2)) goto 90
               k = k + 1
               igsub(k) = kat
               goto 100
 90            kx = kx + 1
               igsubx(kx) = kat
 100        continue
 110     continue
 120  continue

      it = ied2*ncell(1) + 1
      do 160 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 150 ix=ibx,iex
            iix = mod(ix+ncell(1),ncell(1))
            icell = icell0 + iix
            iate = igridadr(icell+1)-1
            do 140 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).gt.gmax(2)) goto 130
               k = k + 1
               igsub(k) = kat
               goto 140
 130           kx = kx + 1
               igsubx(kx) = kat
 140        continue
 150     continue
 160  continue

C     x planes:

      it = ibg1 + 1
      do 200 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 190 iy=iby,iey
            iiy = mod(iy+ncell(2),ncell(2))
            icell = icell0 + iiy*ncell(1)
            iate = igridadr(icell+1)-1
            do 180 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).lt.gmin(1)) goto 170
               k = k + 1
               igsub(k) = kat
               goto 180
 170           kx = kx + 1
               igsubx(kx) = kat
 180        continue
 190     continue
 200  continue

      it = ied1 + 1
      do 240 iz=ibz,iez
         iiz = mod(iz+ncell(3),ncell(3))
         icell0 = iiz*ncell2 + it
         do 230 iy=iby,iey
            iiy = mod(iy+ncell(2),ncell(2))
            icell = icell0 + iiy*ncell(1)
            iate = igridadr(icell+1)-1
            do 220 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).gt.gmax(1)) goto 210
               k = k + 1
               igsub(k) = kat
               goto 220
 210           kx = kx + 1
               igsubx(kx) = kat
 220        continue
 230     continue
 240  continue

      end

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

      subroutine pgetacore(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 getaat(icell, k, igsub)
 10         continue
 20      continue
 30   continue

      end

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

      subroutine pgetaat(icell, k, igsub)

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

      dimension igsub(*)

      iate = igridadr(icell+1)-1
      do 10 iat=igridadr(icell),iate
         k = k + 1
         igsub(k) = igrid(iat)
 10   continue

      end

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

      subroutine pgetainner(i, k, igsub)

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

      dimension igsub(*)


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

      ks = 0

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 getaat(icell, k, igsub)
            enddo
            do iz=izb2,ize2
               iiz = mod(iz+ncell(3),ncell(3))
               icell = icell0 + iiz*ncell2
               call getaat(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 getaat(icell, k, igsub)
            enddo
            do iy=iyb2,iye2
               iiy = mod(iy+ncell(2),ncell(2))
               icell = icell0 + iiy*ncell(1)
               call getaat(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 getaat(icell, k, igsub)
            enddo
            do ix=ixb2,ixe2
              iix = mod(ix+ncell(1),ncell(1))
               icell = icell0 + iix
               call getaat(icell, k, igsub)
            enddo
         enddo
      enddo

      end

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