C
C $Id: gridafnc.F,v 1.3 1998/07/16 16:39:59 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     this file contains all the subroutines needed
C     for atom-wise, grid-based subsetting. These subroutines
C     are called in the driver subroutines atgrid and mixgrid.
C
C     Written by Arjan van der Vaart, July-August 1997
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 getacorners(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:
      logical liyb, lixb, liye, lixe

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

      liyb = ibg(2,i).ge.0
      lixb = ibg(1,i).ge.0
      liye = ied(2,i).lt.ncell(2)
      lixe = ied(1,i).lt.ncell(1)

      if (ibg(3,i).ge.0) then

         if (liyb) then
            icell0 = ibg(3,i)*ncell2 + ibg(2,i)*ncell(1) + 1

            if (lixb) then
               icell = icell0 + ibg(1,i)
               iate = igridadr(icell+1)-1
               do 20 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).lt.gridmin(1,i)) goto 10
                  if (xyz(2,kat).lt.gridmin(2,i)) goto 10
                  if (xyz(3,kat).lt.gridmin(3,i)) goto 10
                  k = k + 1
                  igsub(k) = kat
                  goto 20
 10               kx = kx + 1
                  igsubx(kx) = kat
 20            continue
            endif

            if (lixe) then
               icell = icell0 + ied(1,i)
               iate = igridadr(icell+1)-1
               do 40 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).gt.gridmax(1,i)) goto 30
                  if (xyz(2,kat).lt.gridmin(2,i)) goto 30
                  if (xyz(3,kat).lt.gridmin(3,i)) goto 30
                  k = k + 1
                  igsub(k) = kat
                  goto 40
 30               kx = kx + 1
                  igsubx(kx) = kat
 40            continue
            endif
         endif
         if (liye) then
            icell0 = ibg(3,i)*ncell2 + ied(2,i)*ncell(1) + 1
            if (lixb) then
               icell = icell0 + ibg(1,i)
               iate = igridadr(icell+1)-1
               do 60 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).lt.gridmin(1,i)) goto 50
                  if (xyz(2,kat).gt.gridmax(2,i)) goto 50
                  if (xyz(3,kat).lt.gridmin(3,i)) goto 50
                  k = k + 1
                  igsub(k) = kat
                  goto 60
 50               kx = kx + 1
                  igsubx(kx) = kat
 60            continue
            endif
            if (lixe) then
               icell = icell0 + ied(1,i)
               iate = igridadr(icell+1)-1
               do 80 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).gt.gridmax(1,i)) goto 70
                  if (xyz(2,kat).gt.gridmax(2,i)) goto 70
                  if (xyz(3,kat).lt.gridmin(3,i)) goto 70
                  k = k + 1
                  igsub(k) = kat
                  goto 80
 70               kx = kx + 1
                  igsubx(kx) = kat
 80            continue
            endif
         endif
      endif

      if(ied(3,i).lt.ncell(3)) then
         if (liyb) then
            icell0 = ied(3,i)*ncell2 + ibg(2,i)*ncell(1) + 1
            if (lixb) then
               icell = icell0 + ibg(1,i)
               iate = igridadr(icell+1)-1
               do 100 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).lt.gridmin(1,i)) goto 90
                  if (xyz(2,kat).lt.gridmin(2,i)) goto 90
                  if (xyz(3,kat).gt.gridmax(3,i)) goto 90
                  k = k + 1
                  igsub(k) = kat
                  goto 100
 90               kx = kx + 1
                  igsubx(kx) = kat
 100           continue
            endif
            if (lixe) then
               icell = icell0 + ied(1,i)
               iate = igridadr(icell+1)-1
               do 120 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).gt.gridmax(1,i)) goto 110
                  if (xyz(2,kat).lt.gridmin(2,i)) goto 110
                  if (xyz(3,kat).gt.gridmax(3,i)) goto 110
                  k = k + 1
                  igsub(k) = kat
                  goto 120
 110               kx = kx + 1
                  igsubx(kx) = kat
 120            continue
            endif
         endif
         if (liye) then
            icell0 = ied(3,i)*ncell2 + ied(2,i)*ncell(1) + 1
            if (lixb) then
               icell = icell0 + ibg(1,i)
               iate = igridadr(icell+1)-1
               do 140 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).lt.gridmin(1,i)) goto 130
                  if (xyz(2,kat).gt.gridmax(2,i)) goto 130
                  if (xyz(3,kat).gt.gridmax(3,i)) goto 130
                  k = k + 1
                  igsub(k) = kat
                  goto 140
 130              kx = kx + 1
                  igsubx(kx) = kat
 140           continue
            endif
            if (lixe) then
               icell = icell0 + ied(1,i)
               iate = igridadr(icell+1)-1
               do 160 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).gt.gridmax(1,i)) goto 170
                  if (xyz(2,kat).gt.gridmax(2,i)) goto 170
                  if (xyz(3,kat).gt.gridmax(3,i)) goto 170
                  k = k + 1
                  igsub(k) = kat
                  goto 160
 170              kx = kx + 1
                  igsubx(kx) = kat
 160           continue
            endif
         endif
      endif

      end

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

      subroutine getaribs(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:
      logical lixb, liyb, lizb, lixe, liye, lize

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

      lixb = ibg(1,i).ge.0
      liyb = ibg(2,i).ge.0
      lizb = ibg(3,i).ge.0
      lixe = ied(1,i).lt.ncell(1)
      liye = ied(2,i).lt.ncell(2)
      lize = ied(3,i).lt.ncell(3)

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

C     zx-ribs
      if (lizb.and.lixb) then
         icell0 = ibg(3,i)*ncell2 + ibg(1,i) + 1
         do 30 iy=iyb,iye
            icell = icell0 + iy*ncell(1)
            iate = igridadr(icell+1)-1
            do 20 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).lt.gridmin(1,i)) goto 10
               if (xyz(3,kat).lt.gridmin(3,i)) goto 10
               k = k + 1
               igsub(k) = kat
               goto 20
 10            kx = kx + 1
               igsubx(kx) = kat
 20         continue
 30      continue
      endif

      if (lizb.and.lixe) then
         icell0 = ibg(3,i)*ncell2 + ied(1,i) + 1
         do 60 iy=iyb,iye
            icell = icell0 + iy*ncell(1)
            iate = igridadr(icell+1)-1
            do 50 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).gt.gridmax(1,i)) goto 40
               if (xyz(3,kat).lt.gridmin(3,i)) goto 40
               k = k + 1
               igsub(k) = kat
               goto 50
 40            kx = kx + 1
               igsubx(kx) = kat
 50         continue
 60      continue
      endif

      if (lize.and.lixb) then
         icell0 = ied(3,i)*ncell2 + ibg(1,i) + 1
         do 90 iy=iyb,iye
            icell = icell0 + iy*ncell(1)
            iate = igridadr(icell+1)-1
            do 80 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).lt.gridmin(1,i)) goto 70
               if (xyz(3,kat).gt.gridmax(3,i)) goto 70
               k = k + 1
               igsub(k) = kat
               goto 80
 70            kx = kx + 1
               igsubx(kx) = kat
 80         continue
 90      continue
      endif

      if (lize.and.lixe) then
         icell0 = ied(3,i)*ncell2 + ied(1,i) + 1
         do 120 iy=iyb,iye
            icell = icell0 + iy*ncell(1)
            iate = igridadr(icell+1)-1
            do 110 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).gt.gridmax(1,i)) goto 100
               if (xyz(3,kat).gt.gridmax(3,i)) goto 100
               k = k + 1
               igsub(k) = kat
               goto 110
 100           kx = kx + 1
               igsubx(kx) = kat
 110        continue
 120     continue
      endif

C     zy-ribs

      if (lizb.and.liyb) then
         icell0 = ibg(3,i)*ncell2 + ibg(2,i)*ncell(1) + 1
         do 150 ix=ixb,ixe
            icell = icell0 + ix
            iate = igridadr(icell+1)-1
            do 140 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).lt.gridmin(2,i)) goto 130
               if (xyz(3,kat).lt.gridmin(3,i)) goto 130
               k = k + 1
               igsub(k) = kat
               goto 140
 130           kx = kx + 1
               igsubx(kx) = kat
 140        continue
 150     continue
      endif

      if (lizb.and.liye) then
         icell0 = ibg(3,i)*ncell2 + ied(2,i)*ncell(1) + 1
         do 180 ix=ixb,ixe
            icell = icell0 + ix
            iate = igridadr(icell+1)-1
            do 170 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).gt.gridmax(2,i)) goto 160
               if (xyz(3,kat).lt.gridmin(3,i)) goto 160
               k = k + 1
               igsub(k) = kat
               goto 170
 160           kx = kx + 1
               igsubx(kx) = kat
 170        continue
 180     continue
      endif

      if (lize.and.liyb) then
         icell0 = ied(3,i)*ncell2 + ibg(2,i)*ncell(1) + 1
         do 210 ix=ixb,ixe
            icell = icell0 + ix
            iate = igridadr(icell+1)-1
            do 200 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).lt.gridmin(2,i)) goto 190
               if (xyz(3,kat).gt.gridmax(3,i)) goto 190
               k = k + 1
               igsub(k) = kat
               goto 200
 190           kx = kx + 1
               igsubx(kx) = kat
 200        continue
 210     continue
      endif

      if (lize.and.liye) then
         icell0 = ied(3,i)*ncell2 + ied(2,i)*ncell(1) + 1
         do 240 ix=ixb,ixe
            icell = icell0 + ix
            iate = igridadr(icell+1)-1
            do 230 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(2,kat).gt.gridmax(2,i)) goto 220
               if (xyz(3,kat).gt.gridmax(3,i)) goto 220
               k = k + 1
               igsub(k) = kat
               goto 230
 220           kx = kx + 1
               igsubx(kx) = kat
 230        continue
 240     continue
      endif

C     xy-ribs
      if (liyb.and.lixb) then
         icell0 = ibg(2,i)*ncell(1) + ibg(1,i) + 1
         do 270 iz=izb,ize
            icell = icell0 + iz*ncell2
            iate = igridadr(icell+1)-1
            do 260 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).lt.gridmin(1,i)) goto 250
               if (xyz(2,kat).lt.gridmin(2,i)) goto 250
               k = k + 1
               igsub(k) = kat
               goto 260
 250           kx = kx + 1
               igsubx(kx) = kat
 260        continue
 270     continue
      endif

      if (liyb.and.lixe) then
         icell0 = ibg(2,i)*ncell(1) + ied(1,i) + 1
         do 300 iz=izb,ize
            icell = icell0 + iz*ncell2
            iate = igridadr(icell+1)-1
            do 290 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).gt.gridmax(1,i)) goto 280
               if (xyz(2,kat).lt.gridmin(2,i)) goto 280
               k = k + 1
               igsub(k) = kat
               goto 290
 280           kx = kx + 1
               igsubx(kx) = kat
 290        continue
 300     continue
      endif

      if (liye.and.lixb) then
         icell0 = ied(2,i)*ncell(1) + ibg(1,i) + 1
         do 330 iz=izb,ize
            icell = icell0 + iz*ncell2
            iate = igridadr(icell+1)-1
            do 320 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).lt.gridmin(1,i)) goto 310
               if (xyz(2,kat).gt.gridmax(2,i)) goto 310
               k = k + 1
               igsub(k) = kat
               goto 320
 310           kx = kx + 1
               igsubx(kx) = kat
 320        continue
 330     continue
      endif

      if (liye.and.lixe) then
         icell0 = ied(2,i)*ncell(1) + ied(1,i) + 1
         do 360 iz=izb,ize
            icell = icell0 + iz*ncell2
            iate = igridadr(icell+1)-1
            do 350 iat=igridadr(icell),iate
               kat = igrid(iat)
               if (xyz(1,kat).gt.gridmax(1,i)) goto 340
               if (xyz(2,kat).gt.gridmax(2,i)) goto 340
               k = k + 1
               igsub(k) = kat
               goto 350
 340           kx = kx + 1
               igsubx(kx) = kat
 350        continue
 360     continue
      endif

      end

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

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

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

      dimension igsub(*), igsubx(*)

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

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

C     z planes:

      if (ibg(3,i).ge.0) then
         it = ibg(3,i)*ncell2 + 1
         do 40 iy=iby,iey
            icell0 = iy*ncell(1) + it
            do 30 ix=ibx,iex
               icell = icell0 + ix
               iate = igridadr(icell+1)-1
               do 20 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(3,kat).lt.gridmin(3,i)) goto 10
                  k = k + 1
                  igsub(k) = kat
                  goto 20
 10               kx = kx + 1
                  igsubx(kx) = kat
 20            continue
 30         continue
 40      continue
      endif

      if (ied(3,i).lt.ncell(3)) then
         it = ied(3,i)*ncell2 + 1
         do 80 iy=iby,iey
            icell0 = iy*ncell(1) + it
            do 70 ix=ibx,iex
               icell = icell0 + ix
               iate = igridadr(icell+1)-1
               do 60 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(3,kat).gt.gridmax(3,i)) goto 50
                  k = k + 1
                  igsub(k) = kat
                  goto 60
 50               kx = kx + 1
                  igsubx(kx) = kat
 60            continue
 70         continue
 80      continue
      endif

C     y planes:

      if (ibg(2,i).ge.0) then
         it = ibg(2,i)*ncell(1) + 1
         do 120 iz=ibz,iez
            icell0 = iz*ncell2 + it
            do 110 ix=ibx,iex
               icell = icell0 + ix
               iate = igridadr(icell+1)-1
               do 100 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(2,kat).lt.gridmin(2,i)) goto 90
                  k = k + 1
                  igsub(k) = kat
                  goto 100
 90               kx = kx + 1
                  igsubx(kx) = kat
 100           continue
 110        continue
 120     continue
      endif

      if (ied(2,i).lt.ncell(2)) then
         it = ied(2,i)*ncell(1) + 1
         do 160 iz=ibz,iez
            icell0 = iz*ncell2 + it
            do 150 ix=ibx,iex
               icell = icell0 + ix
               iate = igridadr(icell+1)-1
               do 140 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(2,kat).gt.gridmax(2,i)) goto 130
                  k = k + 1
                  igsub(k) = kat
                  goto 140
 130              kx = kx + 1
                  igsubx(kx) = kat
 140           continue
 150        continue
 160     continue
      endif

C     x planes:

      if (ibg(1,i).ge.0) then
         it = ibg(1,i) + 1
         do 200 iz=ibz,iez
            icell0 = iz*ncell2 + it
            do 190 iy=iby,iey
               icell = icell0 + iy*ncell(1)
               iate = igridadr(icell+1)-1
               do 180 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).lt.gridmin(1,i)) goto 170
                  k = k + 1
                  igsub(k) = kat
                  goto 180
 170              kx = kx + 1
                  igsubx(kx) = kat
 180           continue
 190        continue
 200     continue
      endif

      if (ied(1,i).lt.ncell(1)) then
         it = ied(1,i) + 1
         do 240 iz=ibz,iez
            icell0 = iz*ncell2 + it
            do 230 iy=iby,iey
               icell = icell0 + iy*ncell(1)
               iate = igridadr(icell+1)-1
               do 220 iat=igridadr(icell),iate
                  kat = igrid(iat)
                  if (xyz(1,kat).gt.gridmax(1,i)) goto 210
                  k = k + 1
                  igsub(k) = kat
                  goto 220
 210              kx = kx + 1
                  igsubx(kx) = kat
 220           continue
 230        continue
 240     continue
      endif

      end

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

      subroutine getacore(k, igsub)

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

      dimension igsub(*)

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

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

      do 30 iz=izb,ize
         izn = iz*ncell2 + 1
         do 20 iy=iyb,iye
            icell0 = izn + iy*ncell(1)
            do 10 ix=ixb,ixe
               icell = icell0 + ix
               call getaat(icell, k, igsub)
 10         continue
 20      continue
 30   continue

      end

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

      subroutine getaat(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 getainner(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
      ixb = max(0,ixb)
      ixe = ied(1,i)-1
      ixe = min(ncellm1(1),ixe)
      iyb = ibg(2,i)+1
      iyb = max(0,iyb)
      iye = ied(2,i)-1
      iye = min(ncellm1(2),iye)
      izb1 = ibg(3,i)+1
      izb1 = max(0,izb1)
      ize1 = ibg(3,i-1)-1
      ize1 = min(ncellm1(3),ize1)
      izb2 = ied(3,i-1)+1
      izb2 = max(0,izb2)
      ize2 = ied(3,i)-1
      ize2 = min(ncellm1(3),ize2)

      do iy=iyb,iye
         iyn = iy*ncell(1) + 1
         do ix=ixb,ixe
            icell0 = iyn + ix
            do iz=izb1,ize1
               icell = icell0 + iz*ncell2
               call getaat(icell, k, igsub)
            enddo
            do iz=izb2,ize2
               icell = icell0 + iz*ncell2
               call getaat(icell, k, igsub)
            enddo
         enddo
      enddo

C     y layer

      izb = ibg(3,i-1)
      izb = max(0,izb)
      ize = ied(3,i-1)
      ize = min(ncellm1(3),ize)
      iyb1 = ibg(2,i)+1
      iyb1 = max(0,iyb1)
      iye1 = ibg(2,i-1)-1
      iye1 = min(ncellm1(2),iye1)
      iyb2 = ied(2,i-1)+1
      iyb2 = max(0,iyb2)
      iye2 = ied(2,i)-1
      iye2 = min(ncellm1(2),iye2)
      do ix=ixb,ixe
         ixn = ix + 1
         do iz=izb,ize
            icell0 = ixn + iz*ncell2
            do iy=iyb1,iye1
               icell = icell0 + iy*ncell(1)
               call getaat(icell, k, igsub)
            enddo
            do iy=iyb2,iye2
               icell = icell0 + iy*ncell(1)
               call getaat(icell, k, igsub)
            enddo
         enddo
      enddo

C     x layer


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

      do iy=iyb,iye
         iyn = iy*ncell(1) + 1
         do iz=izb,ize
            icell0 = iyn + iz*ncell2
            do ix=ixb1,ixe1
               icell = icell0 + ix
               call getaat(icell, k, igsub)
            enddo
            do ix=ixb2,ixe2
               icell = icell0 + ix
               call getaat(icell, k, igsub)
            enddo
         enddo
      enddo

      end

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