      SUBROUTINE MAPS

      IMPLICIT CHARACTER ( A-Z )

      COMMON/LSTDAT/LIST,HEAD,MAP
      COMMON/LSTDT2/SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ,MX,MY,MZ

C    *******************************************************************
C    ** ROUTINE TO SET UP A LIST OF NEIGHBOURING CELLS                **
C    *******************************************************************

      INTEGER     NMAX, MMX, MMY, MMZ, NCELL, MAPSIZ
      PARAMETER ( NMAX = 60000 )
      PARAMETER ( MMX=50, MMY=50, MMZ=50, NCELL=MMX*MMY*MMZ )
      PARAMETER ( MAPSIZ = 13*NCELL )

      REAL        SFX,SFY,SFZ,CELLIX,CELLIY,CELLIZ

      INTEGER     LIST(NMAX), HEAD(NCELL), MAP(MAPSIZ), MX, MY,MZ
      INTEGER     IX, IY, IZ, IMAP, ICELL

C    ** STATEMENT FUNCTION TO GIVE CELL INDEX **
      ICELL(IX,IY,IZ) = 1+MOD(IX-1+MX,MX)+MOD(IY-1+MY,MY)*MX
     &                   +MOD(IZ-1+MZ,MZ)*MX*MY

C    ** FIND HALF THE NEAREST NEIGHBOURS OF EACH CELL **
      DO 40 IX = 1, MX
        DO 40 IY = 1, MY
          DO 40 IZ = 1, MZ
            IMAP = (ICELL(IX,IY,IZ)-1)*13
            MAP(IMAP+1 ) = ICELL(IX+1,IY  ,IZ)
            MAP(IMAP+2 ) = ICELL(IX+1,IY+1,IZ)
            MAP(IMAP+3 ) = ICELL(IX  ,IY+1,IZ)
            MAP(IMAP+4 ) = ICELL(IX-1,IY+1,IZ)
            MAP(IMAP+5 ) = ICELL(IX+1,IY  ,IZ-1)
            MAP(IMAP+6 ) = ICELL(IX+1,IY+1,IZ-1)
            MAP(IMAP+7 ) = ICELL(IX  ,IY+1,IZ-1)
            MAP(IMAP+8 ) = ICELL(IX-1,IY+1,IZ-1)
            MAP(IMAP+9 ) = ICELL(IX+1,IY  ,IZ+1)
            MAP(IMAP+10) = ICELL(IX+1,IY+1,IZ+1)
            MAP(IMAP+11) = ICELL(IX  ,IY+1,IZ+1)
            MAP(IMAP+12) = ICELL(IX-1,IY+1,IZ+1)
            MAP(IMAP+13) = ICELL(IX  ,IY  ,IZ+1)
40          CONTINUE

      RETURN
      END
