subroutine geom_vec

! Folgende Groessen muessen gesetzt sein:             B Bunk 7/2010
!     Dimension     ndim                              rev    6/2013
!     Gittergroesse lsize(k), k=1..ndim

! Angelegt und berechnet wird
!     Volumen                    nvol

!     Volumen einer Zelle        nvcell
!     NN-Schritt in der Zelle    nnstep(ib,k), ib=1..nvcell, k=+/-1,..,+/-ndim
!     Permutationsflags dazu     nnflag(ib,k), ib=1..nvcell, k=+/-1,..,+/-ndim

!     Vektorlaenge               lvec
!     NN-Indexfeld im Vektor     nn(i,k), i=1..nvol, k=+/-1,..,+/-ndim

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

! nnstep(ib,k) gibt den Nachbarn des Basispunkts ib bei einem Schritt in
!     Richtung k
!     nnstep(ib,0) ist reserviert.
! nnflag(ib,k) gibt die zugehoerige Permutationsnummer:
!     nnflag(ib,k) = k, wenn der Schritt im Basisblock ueber den Rand springt;
!     nnflag(ib,k) = 0 sonst.
!     nnflag(ib,0) = 0,1 gibt die Farbe von ib (Schachbrettmuster).

! nn(iv,k) gibt, ausgehend vom Vektorindex iv, den Index im Nachbarvektor
!     bei einer Verschiebung (des gesamten Vektors) in Richtung k.
!     nn(iv,0) = iv (keine Permutation).
 
use global

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

do k=1,ndim
   if (mod(lsize(k), 2) /= 0) then
      print '(a,i1,a,i4,a)', &
         'error in geom_vec: lsize(', k,') =', lsize(k), ' is odd'
      stop
   endif
   do lsv=10,1,-1
      if (mod(lsize(k), 2*lsv) == 0) then
         lsvec(k) = lsv
         lscell(k) = lsize(k)/lsv
         exit
      endif
   enddo
enddo
nvol = product(lsize(1:ndim))
nvcell = product(lscell(1:ndim))
lvec = product(lsvec(1:ndim))

if (allocated(nn)) deallocate(nn, nnstep, nnflag)
allocate(nnstep(nvcell, -ndim:ndim))
allocate(nnflag(nvcell, -ndim:ndim))
allocate(nn(lvec, -ndim:ndim))

! Geometriefelder fuer Basiszelle

ibase(1) = 1
do k=1,ndim
   ibase(k+1) = ibase(k)*lscell(k)
enddo
ix(1:ndim) = 1          ! Koord. des Anfangspunkts

do i=1,nvcell           ! Schleife ueber Basispunkte
   nnflag(i,0) = mod(sum(ix(1:ndim)) - ndim, 2) ! Kolorierung
   do k=1,ndim
      nnstep(i,k) = i + ibase(k)                ! Nachbar +k
      nnflag(i,k) = 0
      if (ix(k) == lscell(k)) then
         nnstep(i,k) = nnstep(i,k) - ibase(k+1)
         nnflag(i,k) = k
      endif

      nnstep(i,-k) = i - ibase(k)               ! Nachbar -k
      nnflag(i,-k) = 0
      if (ix(k) == 1) then
         nnstep(i,-k) = nnstep(i,-k) + ibase(k+1)
         nnflag(i,-k) = -k
      endif
   enddo

   do k=1,ndim          ! Koord. des naechsten Punkts
      ix(k) = ix(k) + 1
      if (ix(k) <= lscell(k)) exit
      ix(k) = 1
   enddo
enddo                   ! Ende der Schleife ueber Basispunkte

! Indexfeld fuer Permutation des Vektors

ibase(1) = 1
do k=1,ndim
   ibase(k+1) = ibase(k)*lsvec(k)
enddo
ix(1:ndim) = 1          ! Koord. des Anfangspunkts

do i=1,lvec             ! Schleife ueber Vektorindizes
   nn(i,0) = i                      ! keine Permutation

   do k=1,ndim
      nn(i,k) = i + ibase(k)        ! Nachbar +k
      if (ix(k) == lsvec(k)) nn(i,k) = nn(i,k) - ibase(k+1)

      nn(i,-k) = i - ibase(k)       ! Nachbar -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) <= lsvec(k)) exit
      ix(k) = 1
   enddo
enddo                   ! Ende der Schleife ueber Vektorindizes
end
