SUBROUTINE HMDOT(J,IMERGE,M1,KW,MC,DMS,RNEW,ITERM)
*
*
*       Mass loss from inner hierarchical binary.
*       -----------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
      REAL*8  BODYI(2),W(2)
      REAL*8  M1,MC
      LOGICAL IQUIT
*
*
*       Define merger termination index and relevant KS/c.m. indices.
      ITERM = 0
      IPAIR = KSPAIR
      I1 = 2*IPAIR - 1
      ICM = N + IPAIR
*
*       Include quitting conditions to avoid large changes at end-point.
      IQUIT = .FALSE.
*     IF (RADIUS(J)*SU.GT.500.0) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR.GT.15.0.AND.KSTAR(J).GE.4) IQUIT = .TRUE.
*     IF (ABS(M1 - MC).LT.0.1*M1.OR.MC.GT.1.4) IQUIT = .TRUE.
*     IF (BODY0(J)*ZMBAR - M1.GT.0.5*M1) IQUIT = .TRUE.
*       Note that tidal dissipation needs full KS for rectification in MDOT.
*     IF (KSTARM(IMERGE).LT.0.OR.BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
      IF (BODY(ICM).EQ.0.0D0) IQUIT = .TRUE.
*
*       Quit for mis-identification, advanced evolution or double merger.
      IF (J.LE.0.OR.IQUIT) THEN
          ITERM = 1
          GO TO 50
      END IF
*
*       Update radius for #J and copy #I1 for MDOT & HCORR in case of ghost.
      RSTAR = RNEW
      IF (J.GT.I1) THEN
          RADIUS(J) = RNEW
          RSTAR = RADIUS(I1)
*       Return RNEW for copying back to RADIUS(I1) on ghost & DM/M < 0.01.
      END IF
*
*       Skip further modifications on zero mass loss of same type.
      IF (ABS(DMS/M1).EQ.0.0D0.AND.KW.EQ.KSTAR(J)) GO TO 50
*
*       Set two-body elements for outer orbit and inner semi-major axis.
      SEMI = -0.5*BODY(ICM)/H(IPAIR)
      ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(ICM)*SEMI)
      ECC1 = SQRT(ECC2)
      ZMB0 = CM(1,IMERGE) + CM(2,IMERGE)
      SEMI0 = -0.5*ZMB0/HM(IMERGE)
*
*       Form modified mass ratio and semi-major axes from M*A = const.
      DM = DMS/ZMBAR
      ZMB = ZMB0 - DM
      SEMI1 = SEMI*(ZMB0 + BODY(2*IPAIR))/(ZMB + BODY(2*IPAIR))
      SEMI2 = SEMI0*ZMB0/ZMB
      PMIN = SEMI1*(1.0 - ECC1)
*
*       Evaluate old separation, square regularized velocity, t'' & ECC.
      RI = 0.0D0
      V20 = 0.0
      TD2 = 0.0
      DO 10 K = 1,4
          RI = RI + UM(K,IMERGE)**2
          V20 = V20 + UMDOT(K,IMERGE)**2
          TD2 = TD2 + 2.0*UM(K,IMERGE)*UMDOT(K,IMERGE)
   10 CONTINUE
      ECC2 = (1.0 - RI/SEMI0)**2 + TD2**2/(SEMI0*ZMB0)
      ECC = SQRT(ECC2)
*
*       Determine inclination (use #I1 as first KS component).
      CALL HIMAX(I1,IMERGE,ECC,SEMI0,EMAX,EMIN,ZI,TG,EDAV)
*
*       Evaluate the general stability function (Mardling 2008).
      IF (ECC1.LT.1.0) THEN
          EOUT = ECC1
*       Increase tolerance near sensitive stability boundary (RM 10/2008).
          IF (EOUT.GT.0.9) THEN
              DE = 0.5*(1.0 - EOUT)
              DE = MIN(DE,0.01D0)
*       Add extra amount 0.011 to avoid switching.
              IF (ECC1.GT.0.9) DE = DE + 0.011
              EOUT = EOUT - DE
              PMIN = SEMI1*(1.0 - EOUT)
          END IF
          NST = NSTAB(SEMI2,SEMI1,ECC,EOUT,ZI,CM(1,IMERGE),
     &                                     CM(2,IMERGE),BODY(2*IPAIR))
          IF (NST.EQ.0) THEN
              PCRIT = 0.99*PMIN
              PCR = stability(CM(1,IMERGE),CM(2,IMERGE),BODY(2*IPAIR),
     &                                          ECC,ECC1,ZI)*SEMI2
          ELSE
              PCRIT = 1.01*PMIN
          END IF
      ELSE
          PCRIT = stability(CM(1,IMERGE),CM(2,IMERGE),BODY(2*IPAIR),
     &                                          ECC,ECC1,ZI)*SEMI2
      END IF
*
*       Update pericentre distance on successful stability test or exit.
      IF (PMIN.GT.PCRIT) THEN
          R0(IPAIR) = PCRIT
      ELSE
          ITERM = 1
          GO TO 50
      END IF
*
*       Check for possible circularization (tidal or sequential).
      RP = SEMI2*(1.0 - ECC)
      IF (KZ(27).GT.0.AND.RP.LT.4.0*RADIUS(J).AND.ECC.GT.0.002.AND.
     &    KSTARM(IMERGE).GE.0) THEN
          TC = 1.0D+10
          IF (KZ(27).EQ.2) THEN
              I2 = 0
*       Identify the ghost by searching single bodies.
              DO 12 K = IFIRST,N
                  IF (BODY(K).EQ.0.0D0.AND.
     &                NAME(K).EQ.NAMEG(IMERGE)) THEN
                      I2 = K
                  END IF
   12         CONTINUE
              IF (I2.GT.0) THEN
                  DO 14 K = 1,2
                      BODYI(K) = CM(K,IMERGE)
   14             CONTINUE
                  TG = 1.0
*       Evaluate the circularization time.
                  CALL HICIRC(RP,ECC,I1,I2,BODYI,TG,TC,EC,ED,W)
                  IF (TC.GT.2000) GO TO 18
              END IF
          END IF
          WRITE (6,15)  NAME(J), KSTAR(J), KSTARM(IMERGE), ECC, DMS,
     &                  RP, RADIUS(J), TC
   15     FORMAT (' HMDOT TERM    NAM K* E DM RP R* TC ',
     &                            I6,2I4,F8.4,F7.3,1P,3E10.2)
          ITERM = 1
          GO TO 50
      ELSE IF (ECC.GT.0.002) THEN
          IF (RP.LT.2.0*RADIUS(J)) THEN
              WRITE (6,15)  NAME(J), KSTAR(J), KSTARM(IMERGE), ECC, DMS,
     &                      RP, RADIUS(J)
              ITERM = 1
              GO TO 50
          END IF
      END IF
*
*       Obtain energy change from M*A = const and H = -M/(2*A) (DCH 8/96).
   18 DH = DM/SEMI0*(1.0 - 0.5*DM/ZMB0)
      HM0 = HM(IMERGE)
      HM(IMERGE) = HM(IMERGE) + DH
*
*       Form KS coordinate & velocity scaling factors (general point is OK).
      SEMI2 = -0.5*ZMB/HM(IMERGE)
      C2 = SQRT(SEMI2/SEMI0)
      V2 = 0.5*(ZMB + HM(IMERGE)*RI*(SEMI2/SEMI0))
      C1 = SQRT(V2/V20)
*
*       Re-scale KS variables to new energy (H < 0: constant eccentricity).
      DO 20 K = 1,4
          UM(K,IMERGE) = C2*UM(K,IMERGE)
          UMDOT(K,IMERGE) = C1*UMDOT(K,IMERGE)
*       Note no need to update XREL & VREL for RESTART (c.m. error cancels).
   20 CONTINUE
*
*       Reduce mass of relevant component (ghost is original second member).
      ZMU0 = CM(1,IMERGE)*CM(2,IMERGE)/ZMB0
      KM = 1
      IF (J.GE.IFIRST) KM = 2
      CM(KM,IMERGE) = CM(KM,IMERGE) - DM
*
*       Include corrections to EMERGE & EMDOT (consistency; no net effect!).
      ZMU = CM(1,IMERGE)*CM(2,IMERGE)/ZMB
      DECORR = ZMU*HM(IMERGE) - ZMU0*HM0
      EMERGE = EMERGE + DECORR
      EMDOT = EMDOT - DECORR
*       Compensate for EMERGE contribution in DEGRAV (cf. EVENTS).
      EGRAV = EGRAV - DECORR
*
*       Correct outer orbit for mass loss (use #I1 in case #J is a ghost).
      CALL HCORR(I1,DM,RSTAR)
*
*       Print some diagnostics on non-zero mass loss.
      IF (DMS.GT.1.0D-03) THEN
          WRITE (6,30)  NAME(J), KSTAR(J), IMERGE, KM, M1, DMS, PMIN,
     &                  PCRIT, SEMI2, SEMI1, DECORR
   30     FORMAT (' HMDOT    NAM K* IM KM M1 DM PM PC A A1 DE ',
     &                       I6,3I4,F6.2,1P,5E10.2,0P,F10.6)
      END IF
*
   50 RETURN
*
      END