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