nbsort.f

Go to the documentation of this file.

      SUBROUTINE nbsort(NBL,IBL,NNB,NBLIST)
*
*
*       Sorting of neighbour lists.
*       ---------------------------
*
      include 'common6.h'
      INTEGER  ibl(lmax),nblist(nmax),lp(nmax)
*
*
*       Initialize counter and list pointers.
      nnb = 0
      DO 1 l = 1,nbl
          lp(l) = 2
    1 CONTINUE
*
*       Reset terminator count and minimum particle index.
    5 iend = 0
      jmin = ntot + 1
*
*       Compare current location of all neighbour lists.
      DO 10 ll = 1,nbl
          i = ibl(ll)
          l = lp(ll)
*       Check each list member or count completed searches.
          IF (l.LT.list(1,i) + 2) THEN
              j = list(l,i)
              IF (j.LT.jmin) THEN
                  jmin = j
                  lmin = ll
*       Update pointer for all J <= JMIN (except latest JMIN).
              ELSE IF (j.LE.jmin) THEN
                  lp(ll) = lp(ll) + 1
              END IF
          ELSE
              iend = iend + 1
          END IF
   10 CONTINUE
*
*       Increase pointer & membership and save new index in NBLIST.
      IF (iend.LT.nbl) THEN
          lp(lmin) = lp(lmin) + 1
          nnb = nnb + 1
          nblist(nnb) = jmin
          go to 5
      END IF
*
*       Exit on rare case of one member with zero neighbours (10/2008).
      IF (nnb.EQ.0) go to 70
*
*       Add any integration particles missing from joint neighbour lists.
      new = 0
      fac = nnb/float(ntot)
      l = 1
   20 i = ibl(l)
*       Skip search if body #I is outside the list range.
      IF (nblist(1).GT.i.OR.nblist(nnb).LT.i) go to 50
      IF (nblist(1).EQ.i.OR.nblist(nnb).EQ.i) go to 60
*
*       Search sequential list using incremental expectation values (> 0).
      ig = i*fac
      ig = max(ig,1)
   30 IF (nblist(ig).GT.i) THEN
          inc = (nblist(ig) - i)*fac
          inc = min(inc,3)
          IF (inc.LE.1) THEN
              ig = ig - 1
              IF (nblist(ig).LT.i) go to 50
          ELSE
              ig = ig - inc
              ig = max(ig,1)
              IF (nblist(ig).LT.i) THEN
                  ig = ig + inc - 1
                  IF (nblist(ig).GT.i) THEN
                      ig = ig - 1
                      IF (nblist(ig).GT.i) ig = ig - 1
                  END IF
              ELSE
                  IF (nblist(ig).GT.i) THEN
                      ig = ig - 1
                      IF (nblist(ig).GT.i) ig = ig - 1
                  END IF
              END IF
          END IF
          go to 30
      END IF
*
*       Treat case of NBLIST < I until boundary exceeded or NBLIST > I.
   40 IF (nblist(ig).LT.i) THEN
          inc = (i - nblist(ig))*fac
          inc = max(inc,1)
          ig = ig + inc
          IF (ig.GE.nnb) THEN
              IF (nblist(nnb).EQ.i) go to 60
              go to 50
          END IF
          IF (inc.EQ.1) THEN
              IF (nblist(ig).GT.i) go to 50
          END IF
*       Reduce index by 1 to avoid resonance oscillations (NBLIST > I).
   45     IF (nblist(ig).GT.i) THEN
              ig = ig - 1
              IF (nblist(ig).LT.i) go to 50
              go to 45
          END IF
          go to 40
      END IF
*
*       Skip addition if body #I is already a member (NBLIST(IG) = I).
      go to 60
*
*       Increase counter and add body #I at end of array.
   50 new = new + 1
      nblist(nnb+new) = i
*
*       Continue until last body has been checked.
   60 l = l + 1
      IF (l.LE.nbl) go to 20
*
*       Specify the total membership (Note: new members not sequential).
      nnb = nnb + new
*
   70 RETURN
*
      END