subroutine geom_pbc

! Folgende Groessen muessen gesetzt sein:             B Bunk 11/2002
!     Dimension     ndim                              rev     4/2013
!     Gittergroesse lsize(k), k=1..ndim

! Angelegt und berechnet wird
!     Volumen       nvol
!     NN-Indexfeld  nn(i,k), i=1..nvol, k=+/-1,..,+/-ndim

! nn(i,k) gibt den Index des Nachbarn von i in Richtung k
!     unter Beruecksichtigung periodischer Randbedingungen

use global

integer :: ibase(ndim+1), ix(ndim)

ibase(1) = 1
do k=1,ndim
      ibase(k+1) = ibase(k)*lsize(k)      ! Basis fuer Punktindizes
enddo
nvol = ibase(ndim+1)                      ! Volumen

if (allocated(nn)) deallocate(nn)
allocate(nn(nvol, -ndim:ndim))

ix(1:ndim) = 1          ! Koord. des Anfangspunkts

do i=1,nvol             ! Schleife ueber Punkte
      do k=1,ndim
            nn(i, k) = i + ibase(k)       ! Nachbar x + e_k
            if (ix(k) == lsize(k)) nn(i,k) = nn(i,k) - ibase(k+1)

            nn(i,-k) = i - ibase(k)       ! Nachbar x - e_k
            if (ix(k) == 1) nn(i,-k) = nn(i,-k) + ibase(k+1)
      enddo

      do k=1,ndim       ! Koord. des naechsten Punkts
            ix(k) = ix(k) + 1
            if (ix(k) <= lsize(k)) exit
            ix(k) = 1
      enddo
enddo
end
