SUBROUTINE CHECKL(I,NNB,XI,XIDOT,RS2,FIRR,FREG,FD,FDR) * * * Neighbour list modification. * ---------------------------- * INCLUDE 'common6.h' REAL*8 XI(3),XIDOT(3),DX(3),DV(3),FIRR(3),FREG(3),FD(3),FDR(3) * * ICM = 1 * Only consider high-velocity particles if KZ(37) = 1. IF (KZ(37).EQ.1) GO TO 350 * Omit special treatments of c.m. particles because of force errors. IF (I.GT.N) GO TO 350 * * See whether any neighbour has encounter with body outside sphere. NNB1 = NNB + 1 LJ = 0 * First perform fast test since actual cases are rare. DO 306 L = 2,NNB1 J = ILIST(L) IF (STEP(J).LT.SMIN) LJ = L 306 CONTINUE IF (LJ.EQ.0) GO TO 320 * JCL = ILIST(LJ) ! switched from JCOMP to local variable 4/14. K = 1 308 IF (K.GT.LIST(1,JCL)) GO TO 320 K = K + 1 J = LIST(K,JCL) IF (STEP(J).GT.SMIN) GO TO 308 * Skip single regularized particles or body #I itself. IF (J.LT.IFIRST.OR.J.EQ.I) GO TO 308 A1 = X(1,J) - X(1,JCL) A2 = X(2,J) - X(2,JCL) A3 = X(3,J) - X(3,JCL) RIJ2 = A1**2 + A2**2 + A3**2 * Only accept body #J as neighbour if distance to JCL is < 2*RMIN. IF (RIJ2.GT.RMIN22) GO TO 308 * IF (J.LE.N) GO TO 309 * Also accept pairs satisfying the c.m. force approximation. IF (CMSEP2*R(J-N)**2.GT.RS2) GO TO 308 * * Now see whether body #J has already been included as neighbour. 309 DO 310 L = 1,NNB IF (J.EQ.ILIST(L+1)) GO TO 308 310 CONTINUE * Include body #J in sequential location of neighbour list. L = NNB + 1 312 IF (J.GT.ILIST(L)) GO TO 314 ILIST(L+1) = ILIST(L) * Move other members down by one until appropriate location is free. L = L - 1 IF (L.GT.1) GO TO 312 * 314 ILIST(L+1) = J NNB = NNB + 1 * Only one close encounter neighbour is allowed without test of NNB. NLSMIN = NLSMIN + 1 * * Finally correct irregular and regular force components. DR2 = 0.0 DRDV = 0.0 DO 315 K = 1,3 DX(K) = X(K,J) - XI(K) DV(K) = XDOT(K,J) - XIDOT(K) DR2 = DR2 + DX(K)**2 DRDV = DRDV + DX(K)*DV(K) 315 CONTINUE * DR2I = 1.0/DR2 DR3I = BODY(J)*DR2I*SQRT(DR2I) DRDV = 3.0*DRDV*DR2I * DO 316 K = 1,3 FIRR(K) = FIRR(K) + DX(K)*DR3I FREG(K) = FREG(K) - DX(K)*DR3I FD(K) = FD(K) + (DV(K) - DX(K)*DRDV)*DR3I FDR(K) = FDR(K) - (DV(K) - DX(K)*DRDV)*DR3I 316 CONTINUE * * Include the other component of two recently disrupted pairs. 320 IF (NNB.GE.NNBMAX) GO TO 330 * L = 0 321 L = L + 2 * Advance list index by two for every identified pair. IF (ILIST(L).GT.IFIRST + 3.OR.L.GT.NNB + 1) GO TO 330 * * Set appropriate pair index for either component. JL = ILIST(L) JPAIR = KVEC(JL) IL1 = ILIST(L+1) IPAIR = KVEC(IL1) * Check whether two consecutive list members belong to same pair. IF (JPAIR.EQ.IPAIR) GO TO 321 * The case of index L referring to the last neighbour is permitted. J = 2*JPAIR IF (J.EQ.ILIST(L)) J = J - 1 * Index of the missing component, subject to neighbour test. IF (J.EQ.I) GO TO 323 JCL = ILIST(L) A1 = X(1,JCL) - X(1,J) A2 = X(2,JCL) - X(2,J) A3 = X(3,JCL) - X(3,J) RIJ2 = A1*A1 + A2*A2 + A3*A3 * Only accept #J as neighbour if distance to JCL is < 2*RMIN. IF (RIJ2.LT.RMIN22) GO TO 324 323 L = L - 1 * Increase search index by one only after unsuccessful test. GO TO 321 * * Correct irregular & regular force components. 324 DR2 = 0.0 DRDV = 0.0 DO 325 K = 1,3 DX(K) = X(K,J) - XI(K) DV(K) = XDOT(K,J) - XIDOT(K) DR2 = DR2 + DX(K)**2 DRDV = DRDV + DX(K)*DV(K) 325 CONTINUE * DR2I = 1.0/DR2 DR3I = BODY(J)*DR2I*SQRT(DR2I) DRDV = 3.0*DRDV*DR2I * DO 326 K = 1,3 FIRR(K) = FIRR(K) + DX(K)*DR3I FREG(K) = FREG(K) - DX(K)*DR3I FD(K) = FD(K) + (DV(K) - DX(K)*DRDV)*DR3I FDR(K) = FDR(K) - (DV(K) - DX(K)*DRDV)*DR3I 326 CONTINUE * LJ = NNB + 1 327 IF (J.GT.ILIST(LJ)) GO TO 328 * Move other members down by one until relevant location is free. ILIST(LJ+1) = ILIST(LJ) LJ = LJ - 1 IF (LJ.GT.1) GO TO 327 * 328 ILIST(LJ+1) = J NNB = NNB + 1 NBDIS = NBDIS + 1 * Note that next list index is increased by two after including #J. IF (NNB.LT.NNBMAX) GO TO 321 * * This part includes the other component of an exchanged pair. 330 IF (LISTR(1).EQ.0.OR.NNB.GE.NNBMAX) GO TO 350 * * Do a fast skip if only one disrupted pair in original location. IF (LISTR(1).EQ.2.AND.LISTR(2).LE.IFIRST + 3) GO TO 350 NNB1 = LISTR(1) LJ = 1 + 0.2*FLOAT(NNB) L = 1 332 ICASE = 0 LG = 0 KTIME = 0 334 KTIME = KTIME + 1 335 IF (L.GT.NNB1) GO TO 350 * L = L + 1 JBODY = LISTR(L) IF (JBODY.LE.IFIRST + 3) GO TO 335 * The two last disrupted pairs have already been considered. 336 LG = LG + LJ * First use large increments to speed up the search. IF (LG.GT.NNB + 1) LG = NNB + 1 IF (ILIST(LG).LT.JBODY.AND.LG.LT.NNB + 1) GO TO 336 * LG = LG + 1 338 LG = LG - 1 IF (ILIST(LG).GT.JBODY.AND.LG.GT.2) GO TO 338 IF (ILIST(LG).EQ.JBODY) ICASE = ICASE + KTIME IF (KTIME.EQ.1) GO TO 334 * * See whether any or both of the components have been identified. IF (ICASE.EQ.0.OR.ICASE.EQ.3) GO TO 332 * J = LISTR(L+1-ICASE) * Index of missing component to be included subject to J = I test. JCL = LISTR(L+ICASE-2) * Arguments for J & JCL are L or L - 1 depending on ICASE. A1 = X(1,JCL) - X(1,J) A2 = X(2,JCL) - X(2,J) A3 = X(3,JCL) - X(3,J) RIJ2 = A1*A1 + A2*A2 + A3*A3 * Accept body #J only if distance to JCL is < 2*RMIN (skip J = I). IF (RIJ2.GT.RMIN22.OR.J.EQ.I) GO TO 332 * * Carry out force modifications due to addition of neighbour. 342 DR2 = 0.0 DRDV = 0.0 DO 343 K = 1,3 DX(K) = X(K,J) - XI(K) DV(K) = XDOT(K,J) - XIDOT(K) DR2 = DR2 + DX(K)**2 DRDV = DRDV + DX(K)*DV(K) 343 CONTINUE * DR2I = 1.0/DR2 DR3I = BODY(J)*DR2I*SQRT(DR2I) DRDV = 3.0*DRDV*DR2I * DO 344 K = 1,3 FIRR(K) = FIRR(K) + DX(K)*DR3I FREG(K) = FREG(K) - DX(K)*DR3I FD(K) = FD(K) + (DV(K) - DX(K)*DRDV)*DR3I FDR(K) = FDR(K) - (DV(K) - DX(K)*DRDV)*DR3I 344 CONTINUE * * Include body #J in neighbour list and increase NNB. 345 K = NNB + 1 346 IF (J.GT.ILIST(K)) GO TO 348 * Move other members down by one until relevant location is free. ILIST(K+1) = ILIST(K) K = K - 1 IF (K.GT.1) GO TO 346 * 348 ILIST(K+1) = J NNB = NNB + 1 NBDIS2 = NBDIS2 + 1 * * Check list of high velocity intruders for early retention. 350 IF (LISTV(1).EQ.0.OR.ICM.EQ.-1) GO TO 355 L = 2 351 J = LISTV(L) IF (J.EQ.I) GO TO 354 A1 = X(1,J) - XI(1) A2 = X(2,J) - XI(2) A3 = X(3,J) - XI(3) RIJ2 = A1**2 + A2**2 + A3**2 * IF (RIJ2.GT.4.0*RS2.OR.RIJ2.LT.RS2) GO TO 354 * Simplify to include GPU (standard code needs a few extra searches). IF (RIJ2.GT.4.0*RS2) GO TO 354 A4 = XDOT(1,J) - XDOT(1,I) A5 = XDOT(2,J) - XDOT(2,I) A6 = XDOT(3,J) - XDOT(3,I) A7 = A1*A4 + A2*A5 + A3*A6 IF (A7.GT.0.0) GO TO 354 P2 = RIJ2 - A7**2/(A4**2 + A5**2 + A6**2) * Accept if impact parameter < RS. IF (P2.GT.RS2) GO TO 354 * See if body #J has been included by other procedures. DO 353 K = 1,NNB IF (J.EQ.ILIST(K+1)) GO TO 354 353 CONTINUE ICM = -1 * Redundant indicator (normally > 0) used for joint procedure. NBDIS2 = NBDIS2 - 1 NBFAST = NBFAST + 1 * * Distinguish betweeen single particle treatment and perturbed case. IF (NNB.LE.NNBMAX) THEN IF (I.LE.N) GO TO 342 RIJ2 = (XI(1) - X(1,J))**2 + (XI(2) - X(2,J))**2 + & (XI(3) - X(3,J))**2 * Adopt c.m. approximation if body #J is also single. IF (RIJ2.GT.CMSEP2*R(I-N)**2.AND.J.LE.N) GO TO 342 ELSE GO TO 355 END IF * * Set KS component indices (also body #J if inside c.m. approximation). I1 = 2*(I - N) - 1 I2 = I1 + 1 J1 = J J2 = 0 IF (J.GT.N) THEN IF (RIJ2.LT.CMSEP2*R(J-N)**2) THEN J1 = 2*(J - N) - 1 J2 = J1 + 1 END IF END IF * * Begin evaluation of all relevant interactions (at most 4 terms). K = J1 360 L = I1 362 A1 = X(1,K) - X(1,L) A2 = X(2,K) - X(2,L) A3 = X(3,K) - X(3,L) RIJ2 = A1*A1 + A2*A2 + A3*A3 DV(1) = XDOT(1,K) - XDOT(1,L) DV(2) = XDOT(2,K) - XDOT(2,L) DV(3) = XDOT(3,K) - XDOT(3,L) * * Employ the appropriate mass-weighted expression. DR2I = 1.0/RIJ2 DR3I = BODY(K)*BODY(L)*DR2I*SQRT(DR2I)/BODY(I) DRDV = 3.0*(A1*DV(1) + A2*DV(2) + A3*DV(3))*DR2I * * Add irregular F & FDOT and subtract regular terms. FIRR(1) = FIRR(1) + A1*DR3I FIRR(2) = FIRR(2) + A2*DR3I FIRR(3) = FIRR(3) + A3*DR3I FD(1) = FD(1) + (DV(1) - A1*DRDV)*DR3I FD(2) = FD(2) + (DV(2) - A2*DRDV)*DR3I FD(3) = FD(3) + (DV(3) - A3*DRDV)*DR3I FREG(1) = FREG(1) - A1*DR3I FREG(2) = FREG(2) - A2*DR3I FREG(3) = FREG(3) - A3*DR3I FDR(1) = FDR(1) - (DV(1) - A1*DRDV)*DR3I FDR(2) = FDR(2) - (DV(2) - A2*DRDV)*DR3I FDR(3) = FDR(3) - (DV(3) - A3*DRDV)*DR3I * * Consider each interaction in turn (body #J may be resolved). L = L + 1 IF (L.EQ.I2) GO TO 362 K = K + 1 IF (K.EQ.J2) GO TO 360 * * Include body #J in the neighbour list. GO TO 345 * 354 L = L + 1 IF (L.LE.LISTV(1) + 1) GO TO 351 * 355 RETURN * END