SUBROUTINE SUBINT(IQ,I10)
*
*
*       Decision-making for subsystems.
*       -------------------------------
*
      INCLUDE 'common6.h'
      COMMON/CLUMP/   BODYS(NCMAX,5),T0S(5),TS(5),STEPS(5),RMAXS(5),
     &                NAMES(NCMAX,5),ISYS(5)
      COMMON/IHIST/ IBIN(1000)
      COMMON/KSPAR/ ISTAT(KMAX)
      INTEGER LISTKS(KMAX),NEXTL(KMAX)
      PARAMETER (NLMAX=3)
      REAL*8  TSLIST(KMAX)
      SAVE IB,LY,NK
      DATA IB,LY /0,0/
*     SAVE COMP,TT1
*     DATA COMP /0.0D0/
*
      TBLIST = TIME
*
*     TT1 = DBLTIM()  ! Needs wtime.o from GPU2 (also in Makefile)
*       Initialize the parallel KS counters.
      IF (IB.EQ.0) THEN
          DO 100 K = 1,1000
              IBIN(K) = 0
  100     CONTINUE
          IB = 1
          DTB = 0.0
      END IF
*
*       See whether to advance any KS solutions at start of block-step.
    1 IF (NPAIRS.GT.0) THEN
*     Obtain list of all KS pairs due in interval DTB.
          IF (TBLIST.LE.TBLOCK.OR.NPAIRS.NE.NBPREV) THEN
              IF (DTB.EQ.0.0D0.OR.DTB.GT.1.0D+06) THEN
                  DTB = MAX(DTMIN,TBLOCK - TPREV)
              END IF
    2         TBLIST = TPREV + DTB
              TBLIST = MAX(TBLOCK,TBLIST)
              NNTB = 0
              DO 3 JPAIR = 1,NPAIRS
                  J1 = 2*JPAIR - 1
                  IF (T0(J1) + STEP(J1).LE.TBLIST) THEN
                      NNTB = NNTB + 1
                      KBLIST(NNTB) = J1
                      TSLIST(NNTB) = T0(J1) + STEP(J1)
                  END IF
    3         CONTINUE
*       Increase interval on zero membership.
              IF (NNTB.EQ.0) THEN
                  DTB = 2.0*DTB
                  GO TO 2
              END IF
*       Stabilize interval on membership of 2*SQRT(NPAIRS).
              NBTRY = 2*SQRT(FLOAT(NPAIRS))
              IF (NNTB.GT.NBTRY)  DTB = 0.75*DTB
              IF (NNTB.LT.NBTRY)  DTB = 1.25*DTB
          END IF
*
*       Determine quantized value of TMIN (next small block-step).
    5     TMIN = TBLOCK
          NK = 0
          DO 4 L = 1,NNTB
              I1 = KBLIST(L)
*             IF (TSLIST(I1).LE.TBLOCK) THEN
              IF (TSLIST(L).LE.TBLOCK) THEN
                  NK = NK + 1
                  LISTKS(NK) = I1
              END IF
    4     CONTINUE
          IF (NK.EQ.0) THEN
              TIME = TBLOCK
              GO TO 60
          END IF
*
          DO 6 L = 1,NK
            I1 = LISTKS(L)
*           TMIN = MIN(TSLIST(I1),TMIN)
            TMIN = MIN(TSLIST(L),TMIN)
    6     CONTINUE
*
*       Update time and form list of binaries due at TMIN.
          TIME = TMIN
          LENGTH = 0
          DO 8 L = 1,NK
            I1 = LISTKS(L)
            IF (T0(I1) + STEP(I1).EQ.TMIN) THEN
              LENGTH = LENGTH + 1
              NEXTL(LENGTH) = I1
            END IF
    8     CONTINUE
*
*       Check possible chain integration on zero LENGTH.
          IF (LENGTH.EQ.0) THEN
            GO TO 30 
          END IF
*
*       Clean current length of ISTAT.
          DO 9 L = 1,LENGTH
            ISTAT(L) = 0
    9     CONTINUE
*
*       Begin KS integration (sequential or parallel).
          IF (LENGTH.LE.NLMAX) THEN
            DO 10 LI = 1,LENGTH
              I1 = NEXTL(LI)
              CALL KSINT(I1,LI)
   10       CONTINUE
            NBPRED = NBPRED + LENGTH
*
*       Search non-zero flags (ISTAT < 0 means collision).
            IF (IQ.EQ.0) THEN
              DO 12 LI = 1,LENGTH
                IF (ISTAT(LI).NE.0) THEN
                  IQ = ISTAT(LI)
                  I10 = NEXTL(LI)
                  KSPAIR = KVEC(I10)
*       Note IPHASE must be defined for chain (not needed for collision).
                  IPHASE = IQ
                  GO TO 15
                END IF 
   12         CONTINUE
            END IF
*
*       Treat collision explicitly before quitting (suppressed in KSINT).
   15       IF (IQ.LT.0) THEN
              IP = KVEC(I10)
              KSPAIR = IP
              QPERI = R(IP)
              IQCOLL = -2
              CALL CMBODY(QPERI,2)
*       Ensure block-step time for normal continuation.
              TIME = TBLOCK
*       Note rest of the block will be done next time (need GO TO 999 first).
              GO TO 60
            END IF
*       Complete all block-steps before special procedure (except coll).
            IF (TMIN.LT.TBLOCK.AND.IQ.EQ.0) GO TO 5
*
            TIME = TBLOCK
          ELSE
*
*       Perform the parallel execution loop.
!$omp parallel do private(LI,I1)
            DO 20 LI = 1,LENGTH
              I1 = NEXTL(LI)
              CALL KSINTP(I1,LI)
   20       CONTINUE
!$omp end parallel do
            NSTEPU = NSTEPU + LENGTH
*       Note that unperturbed KS are now included in parallel step counter.
*
      IF (N.GT.0) GO TO 90
            LB = FLOAT(LENGTH)/FLOAT(NLMAX)
            LB = MAX(LB,1)
            LY = MAX(LY,LB)
            LY = MIN(LY,20)
            IBIN(LB) = IBIN(LB) + LENGTH
            IF (ABS(TNEXT - TBLOCK).LT.1.0D-04) THEN
              WRITE (6,22)  (IBIN(K),K=1,LY)
   22         FORMAT (' IBIN  ',20I9)
              ISUM = 0
              DO 23 K = 1,LY
                ISUM = ISUM + IBIN(K)
   23         CONTINUE
              WRITE (6,24)  LY, ISUM
   24         FORMAT (' LY ISUM  ',I4,I10)
            END IF
   90   CONTINUE
*
*       Search non-zero flags (ISTAT < 0 for collision).
            IF (IQ.EQ.0) THEN
              DO 25 LI = 1,LENGTH
                IF (ISTAT(LI).NE.0) THEN
                  IQ = ISTAT(LI)
                  I10 = NEXTL(LI)
                  IPHASE = IQ
                  KSPAIR = KVEC(I10)  ! might be needed somewhere.
                  GO TO 28
                END IF 
   25         CONTINUE
            END IF
*
*       Treat collision explicitly before quitting (suppressed in KSINTP).
   28       IF (IQ.LT.0) THEN
              IP = KVEC(I10)
              KSPAIR = IP
              QPERI = R(IP)
              IQCOLL = -2
              CALL CMBODY(QPERI,2)
*       Note rest of the block will be done next time (need GO TO 999 first).
              GO TO 60
            END IF
*
*       Check more block-step levels unless collision (final exit at TBLOCK).
            IF (TMIN.LT.TBLOCK.AND.IQ.EQ.0) GO TO 5
          END IF
*       Copy original block time at end of KS treatment.
          TIME = TBLOCK
          NBPREV = NPAIRS
      END IF
*
*       Check time for advancing any triple, quad or chain regularization.
   30 IF (NSUB.GT.0) THEN
   40     TSUB = 1.0D+10
          DO 45 L = 1,NSUB
              IF (TS(L).LT.TSUB) THEN
                  ISUB = L
                  TSUB = TS(L)
              END IF
   45     CONTINUE
*
          IF (TSUB.LE.TBLOCK) THEN
              TIME = TSUB
*       Decide between triple, quad or chain.
              IF (ISYS(ISUB).EQ.1) THEN
*       Update unperturbed size of subsystem and copy c.m. step.
                  CALL EXTEND(ISUB)
                  CALL TRIPLE(ISUB)
              ELSE IF (ISYS(ISUB).EQ.2) THEN
                  CALL EXTEND(ISUB)
                  CALL QUAD(ISUB)
              ELSE
                  IF (STEPS(ISUB).LT.0.0D0) THEN
                      STEPS(ISUB) = 1.0D-10
                      GO TO 50
                  END IF
                  CALL CHAIN(ISUB)
                  IF (ISUB.GT.0) THEN
                      IF (STEPS(ISUB).LT.0.0D0) THEN
                          STEPS(ISUB) = 1.0D-10
                          GO TO 50
                      END IF
                  END IF
              END IF
*
*       Check for termination (set TPREV < TIME and set IQ < 0).
              IF (ISUB.LT.0.OR.IPHASE.LT.0) THEN
                  TPREV = TIME - STEP(NTOT)
                  IQ = -1
              END IF
              GO TO 40
          END IF
   50     TIME = TBLOCK
      END IF
*
*     TT2 = DBLTIM()
*     COMP = COMP + (TT2 - TT1)
*     IF (ABS(TIME-TCRIT).LT.1.0D-05) WRITE (6,70)  COMP
*  70 FORMAT (' SUBINT   ',1P,E12.3)
*
   60 RETURN
*
      END