checkl.f

Go to the documentation of this file.

      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
*
      jcomp = ilist(lj)
      k = 1
  308 IF (k.GT.list(1,jcomp)) go to 320
      k = k + 1
      j = list(k,jcomp)
      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,jcomp)
      a2 = x(2,j) - x(2,jcomp)
      a3 = x(3,j) - x(3,jcomp)
      rij2 = a1**2 + a2**2 + a3**2
*       Only accept body #J as neighbour if distance to JCOMP 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
      jcomp = ilist(l)
      a1 = x(1,jcomp) - x(1,j)
      a2 = x(2,jcomp) - x(2,j)
      a3 = x(3,jcomp) - x(3,j)
      rij2 = a1*a1 + a2*a2 + a3*a3
*       Only accept #J as neighbour if distance to JCOMP 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.
      jcomp = listr(l+icase-2)
*       Arguments for J & JCOMP are L or L - 1 depending on ICASE.
      a1 = x(1,jcomp) - x(1,j)
      a2 = x(2,jcomp) - x(2,j)
      a3 = x(3,jcomp) - x(3,j)
      rij2 = a1*a1 + a2*a2 + a3*a3
*       Accept body #J only if distance to JCOMP 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