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