SUBROUTINE SUBSYS(NSYS,CM)
*
*
*       Initialization of subsystem.
*       ----------------------------
*
      INCLUDE 'common6.h'
      COMMON/CLUMP/   BODYS(NCMAX,5),T0S(5),TS(5),STEPS(5),RMAXS(5),
     &                NAMES(NCMAX,5),ISYS(5)
      REAL*8  CM(7)
*
*
*       Increase subsystem counter and set current index.
      NSUB = NSUB + 1
      ISUB = NSUB
*
*       Set zero name & mass in last component to distinguish triple case.
      IF (NSYS.EQ.3) THEN
          NAMES(4,ISUB) = 0
          BODYS(4,ISUB) = 0.0D0
      END IF
*
*       Save global masses & names of new subsystem components.
      BX = 0.0
      DO 10 L = 1,NSYS
          J = JLIST(L)
          BODYS(L,ISUB) = BODY(J)
          NAMES(L,ISUB) = NAME(J)
          IF (BODY(J).GT.BX) THEN
              BX = BODY(J)
              LX = L
          END IF
   10 CONTINUE
*
*       Form ghosts and initialize integration variables (NSYS = 3 or 4).
      DO 20 L = 1,NSYS
          J = JLIST(L)
          BODY(J) = 0.0D0
          T0(J) = 1.0E+06
          LIST(1,J) = 0
          DO 15 K = 1,3
              X0DOT(K,J) = 0.0D0
              XDOT(K,J) = 0.0D0
              F(K,J) = 0.0D0
              FDOT(K,J) = 0.0D0
              D0(K,J) = 0.0D0
              D1(K,J) = 0.0D0
              D2(K,J) = 0.0D0
              D3(K,J) = 0.0D0
              D0R(K,J) = 0.0D0
              D1R(K,J) = 0.0D0
              D2R(K,J) = 0.0D0
              D3R(K,J) = 0.0D0
   15     CONTINUE
*       Set large X0 & X to avoid perturber selection (no escape).
          X0(1,J) = 1.0E+06
          X(1,J) = 1.0E+06
   20 CONTINUE
*
*       Place c.m. of subsystem in first location (ICOMP may switch!).
      ICOMP = JLIST(1)
*       Ensure heaviest body is selected in case of ARchain.
      IF (ISYS(ISUB).EQ.4) THEN
          ICOMP = JLIST(LX)  
      END IF
*       Define zero name for identification (only chain or ARchain c.m.).
      IF (ISYS(ISUB).GE.3) NAME(ICOMP) = 0
*
*       Initialize c.m in ICOMP (redefined above).
      T0(ICOMP) = TIME
      BODY(ICOMP) = CM(7)
      DO 30 K = 1,3
          X(K,ICOMP) = CM(K)
          X0(K,ICOMP) = CM(K)
          XDOT(K,ICOMP) = CM(K+3)
          X0DOT(K,ICOMP) = CM(K+3)
   30 CONTINUE
*
*       Predict coordinates & velocities for all neighbours (order FDOT).
      NNB = LIST(1,ICOMP)
      CALL XVPRED(ICOMP,NNB)
*
*       Obtain new neighbour list (to ensure membership > 0).
      RS0 = RS(ICOMP)
      CALL NBLIST(ICOMP,RS0)
*
*       Construct force polynomials for c.m. motion.
      CALL FPOLY1(ICOMP,ICOMP,0)
      CALL FPOLY2(ICOMP,ICOMP,0)
*
*       Initialize decision-making variables for multiple regularization.
      T0S(ISUB) = TIME
      TS(ISUB) = TIME
      STEPS(ISUB) = STEP(ICOMP)
*
*       Obtain maximum unperturbed separation based on dominant neighbour.
      CALL EXTEND(ISUB)
*
      RETURN
*
      END