roche.f

Go to the documentation of this file.

***
      SUBROUTINE roche(IPAIR)
*
*
*       Roche lobe overflow.
*       Updated 6/1/98 & 28/11/01 by J. Hurley.
*       ---------------------------------------
*
      include 'common6.h'
      common/rche/ cmrch(13,mmax),namer(2,mmax),kstarr(3,mmax)
      REAL*8 tscls(20),lums(10),gb(10),tm,tn,tphys0
      REAL*8 m01,m1,age,r1,lum,mc,rcc,sep,m1ce,m2ce,rce,km0,km,m2,mcx
      REAL*8 mass0(2),mass(2),massc(2),rad(2),radx(2),rol(2),
     &       radc(2),q(2),aj(2),tkh(2),tms(2),tbgb(2),
     &       jspin(2),ospin(2),djspin(2),dtmi(2)
      REAL*8 menv(2),renv(2),k2str(2)
      REAL*8 jorb,oorb,rlperi,djorb,djmb,djgr,djt,rdmin,rdisk
      REAL*8 jspbru,ospbru
      REAL*8 dms(2),dma(2),dmr(2),ivsqm,vwind2,vorb2
      REAL*8 dm1,dm2,dm22,dspin,dspin0,delet,delet1,eqspin,mew,mt2
      REAL*8 tc,tsyn,ttid,fac,fac0,ecc0,dtmmin,djtk(2),dspink(2)
      REAL*8 mch,epsnov,eddfac,gamm1
      parameter(mch=1.44d0,epsnov=0.d0,eddfac=100.d0,gamm1=-1.d0)
      REAL*8 k2,k3,acc1,acc2,beta,xi,aursun
      parameter(k3=0.21d0,acc1=3.920659d+08,acc2=1.5d0)
      parameter(beta=0.125d0,xi=1.d0,aursun=214.95d0)
      REAL*8 mlwind,rl
      EXTERNAL mlwind,rl
      LOGICAL coals,disk,novae,supedd,isave
      CHARACTER*5 ch5
      SAVE iwarn,igr,imb
      DATA iwarn,igr,imb /0,0,0/
      SAVE first
      LOGICAL first
      DATA first /.true./
*
*
*       Set components & c.m. index and initialize counters & interval.
      i1 = 2*ipair - 1
      i2 = i1 + 1
      i = n + ipair
      iter = 0
      iqcoll = 0
      dt1 = 0.0
      isave = .false.
      jkick = 0
      inew = 0
*
*       Form semi-major axis and eccentricity.
      semi = -0.5d0*body(i)/h(ipair)
      ecc2 = (1.d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*body(i))
      ecc = sqrt(ecc2)
*
*       Set age, mass ratio, Roche radius & radius (SU) for each star.
      j = i1
      DO 5 k = 1,2
          q(k) = body(j)/(body(i) - body(j))
          rol(k) = rl(q(k))*semi*su
          rad(k) = radius(j)*su
          dms(k) = 0.d0
          j = i2
    5 CONTINUE
*
*       Determine indices for primary & secondary star (donor & accretor).
      IF (rad(1)/rol(1).GE.rad(2)/rol(2)) THEN
          j1 = i1
          j2 = i2
          rl1 = rol(1)
      ELSE
          j1 = i2
          j2 = i1
          rl1 = rol(2)
          rr = rad(1)
          rad(1) = rad(2)
          rad(2) = rr
          rr = rol(1)
          rol(1) = rol(2)
          rol(2) = rr
          rr = q(1)
          q(1) = q(2)
          q(2) = rr
      END IF
*
*       Exit if physical radius is smaller than Roche radius.
      IF (rad(1).LE.rl1) THEN
          go to 200
      ELSE
*       Set the time to the latest evolution time of the components.
          tphys0 = max(tev0(j1),tev0(j2))*tstar
          tphys = tphys0 + toff*tstar
*       Save the original stellar types and initialize coalescence flag.
          kw1 = kstar(j1)
          kw2 = kstar(j2)
          coals = .false.
*         write(6,*)' start roche ',name(j1),name(j2),kw1,kw2,
*    &                              rad(1)/rl1,rad(2)/rol(2)
*
*       Set separation & masses (solar units) and mass ratios.
          sep = semi*su
          semi0 = semi
          mass(1) = body(j1)*zmbar
          mass(2) = body(j2)*zmbar
          mass0(1) = body0(j1)*zmbar
          mass0(2) = body0(j2)*zmbar
*
*       Set the orbital spin and angular momentum.
          tb = yrs*semi*sqrt(semi/body(i))
          tk = days*tb/yrs
          oorb = twopi/tb
          jorb = mass(1)*mass(2)/(mass(1)+mass(2))
     &           *sqrt(1.d0-ecc*ecc)*sep*sep*oorb
*
          dm1 = 0.d0
          dm2 = 0.d0
*
*       Initialize local arrays.
          j = j1
          DO 97 k = 1,2
             m01 = mass0(k)
             m1 = mass(k)
             mc = 0.d0
             kw = kstar(j)
             age = tev0(j)*tstar - epoch(j)
             CALL star(kw,m01,m1,tm,tn,tscls,lums,gb,zpars)
             CALL hrdiag(m01,age,m1,tm,tn,tscls,lums,gb,zpars,
     &                   r1,lum,kw,mc,rcc,m1ce,rce,k2)
             massc(k) = mc
             aj(k) = age
             menv(k) = m1ce
             renv(k) = rce
             k2str(k) = k2
             radx(k) = r1
             IF(k.EQ.1)THEN
                radx(k) = min(radx(k),rol(k))
                radx(k) = max(rcc,radx(k))
             ENDIF
             jspin(k) = max(spin(j)*spnfac,1.0d-10)
             ospin(k) = jspin(k)/(k2*radx(k)*radx(k)*(m1-mc) + 
     &                  k3*rcc*rcc*mc)
*            IF(K.EQ.1)THEN
*               IF(OSPIN(1).GT.OORB)THEN
*       Force co-rotation of primary and orbit to ensure that the
*       tides do not lead to unstable Roche.
*                  OSPIN(1) = OORB
*                  JSPIN(1) = OSPIN(1)*(K2*R1*R1*(M1-MC)+K3*RCC*RCC*MC)
*               ENDIF
*            ENDIF
             rlperi = rol(k)*(1.d0 - ecc)
             dmr(k) = mlwind(kw,lum,r1,m1,mc,rlperi,zmet)
             dma(k) = 0.d0
             tms(k) = tm
             tbgb(k) = tscls(1)
             j = j2
   97     CONTINUE
*
*       Implement any mass loss between TEV0 and the Roche starting time.
*       This is a check for the case where one star is catching up on the
*       other and Roche occurs in the meantime!
          j = j1
          vorb2 = acc1*(mass(1)+mass(2))/sep
          ivsqm = 1.d0/sqrt(1.d0-ecc*ecc)
          DO 98 k = 1,2
             dt = 1.0d+06*(tphys0 - tev0(j)*tstar)
             IF(dt.EQ.0.0) goto 99
*            VWIND2 = 2.D0*BETA*ACC1*MASS(3-K)/RAD(3-K)
*            OMV2 = (1.D0 + VORB2/VWIND2)**(3.D0/2.D0)
*            DMA(K) = IVSQM*ACC2*DMR(3-K)*((ACC1*MASS(K)/VWIND2)**2)/
*    &                (2.D0*SEP*SEP*OMV2)
*            DMA(K) = MIN(DMA(K),0.8D0*DMR(3-K))
             dms(k) = (dmr(k) - dma(k))*dt
             djspin(k) = (2.d0/3.d0)*(dmr(k)*radx(k)*radx(k)*ospin(k) -
     &                    xi*dma(k)*radx(3-k)*radx(3-k)*ospin(3-k))*dt
             kw = kstar(j)
             IF(kw.LT.10)THEN
                dms(k) = min(dms(k),mass(k)-massc(k))
                CALL magbrk(kw,mass(k),menv(k),rad(k),ospin(k),djmb)
                djspin(k) = djspin(k) + djmb*dt
             ENDIF
             mass(k) = mass(k) - dms(k)
             jspin(k) = max(jspin(k) - djspin(k),1.0d-10)
             IF(kw.LE.2.OR.kw.EQ.7)THEN
                m01 = mass0(k)
                mass0(k) = mass(k)
                body0(j) = mass(k)/zmbar
                CALL star(kw,mass0(k),mass(k),tm,tn,tscls,lums,gb,zpars)
                IF(kw.EQ.2)THEN
                   IF(gb(9).LT.massc(k).OR.m01.GT.zpars(3))THEN
                      mass0(k) = m01
                   ELSE
                      epoch(j) = tm + (tscls(1)-tm)*(aj(k)-tms(k))/
     &                                (tbgb(k)-tms(k))
                      epoch(j) = tev0(j)*tstar - epoch(j)
                   ENDIF
                ELSE
                   epoch(j) = tev0(j)*tstar - aj(k)*tm/tms(k)
                ENDIF
             ELSE
                CALL star(kw,mass0(k),mass(k),tm,tn,tscls,lums,gb,zpars)
             ENDIF
             age = tphys0 - epoch(j)
             IF((age-tn).GT.1.0d-02)THEN
                IF(kw.LE.6)THEN
                   age = min(age,0.9999*tscls(11))
                ELSE
                   age = min(age,0.9999*tscls(8))
                ENDIF
                epoch(j) = tphys0 - age
             ENDIF
   99        j = j2
   98     CONTINUE
*
          sep = sep*(mass(1)+mass(2)+dms(1)+dms(2))/(mass(1)+mass(2))
          semi = sep/su
          tb = yrs*semi*sqrt(semi*zmbar/(mass(1)+mass(2)))
          tk = days*tb/yrs
          oorb = twopi/tb
          jorb = mass(1)*mass(2)/(mass(1)+mass(2))
     &           *sqrt(1.d0-ecc*ecc)*sep*sep*oorb
*
*       Define first/second Roche overflow stage and increase event counter.
          IF (kstar(i).LE.10.OR.
     &        (kstar(i).GT.10.AND.mod(kstar(i),2).EQ.0)) THEN
              kstar(i) = kstar(i) + 1
              nro = nro + 1
              inew = 1
              WRITE (6,8)  name(j1), name(j2), tphys, kw1, kw2,
     &                     kstar(i), mass, sep, tk, rad, rl1
    8         FORMAT (' NEW ROCHE    NAM TP K* M A P R* RL ',
     &                               2i6,f9.2,3i4,2f8.4,f8.2,1p,4e10.2)
              IF(kstar(i).EQ.50)THEN
                 WRITE(6,9)name(j1),name(j2),kw1,kw2
    9            FORMAT(' WARNING: TOO MUCH ROCHE  NAM K* ',2i6,2i3)
              ENDIF
              igr = 0
              imb = 0
              IF (kz(8).GT.3) THEN
                  CALL binev(ipair)
              END IF
*
              IF(first)THEN
                 OPEN(unit=85,status='NEW',form='FORMATTED',
     &                file='ROCHE')
                 first = .false.
                 WRITE (85,94)
   94            FORMAT (/,' NAM1  NAM2  K1 K2   TPHYS     AGE1     ',
     &                     ' AGE2     M01    M02    M1     M2     Z ',
     &                     '     e        P        JSPIN1      JSPIN2')
              ENDIF
              ch5 = ' NEW '
              WRITE(85,95)name(j1),name(j2),kstar(j1),kstar(j2),
     &                    kstar(i),tphys,aj(1),aj(2),
     &                    mass0(1),mass0(2),mass(1),mass(2),
     &                    zmet,ecc,tk,jspin(1),jspin(2),ch5
   95         FORMAT(2i7,3i3,3f10.3,4f7.3,f7.4,f6.3,1p,3e12.4,a5)
              CALL flush(85)
          ENDIF
*
*       Evaluate 1/10 of nuclear time interval (yrs).
          m1 = mass(1)
          CALL trdot(j1,dtm,m1)
          dtm = dtm*tstar
          dtmi(1) = dtm
          dtmi(2) = dtmi(1)
*
*       Include time-scale check for GR & MB braking for small separations.
*       (Note: use old orbital expressions)
          IF (sep.LT.10.0) THEN
              dsepg = -1.663d-09*mass(1)*mass(2)*(mass(1) + mass(2))/
     &                                          (sep*sep*sep)*1.0d+06
              dtgr = -0.1d0*sep/dsepg
              dtm = min(dtgr,dtm)
              IF (kstar(j1).LE.1.AND.
     &            mass(1).GT.0.35.AND.mass(1).LT.1.25) THEN
                  dsepm = -4.574d-07*(mass(1) + mass(2))**2*
     &                        radx(1)**3/(mass(2)*sep**4)*1.0d+06
                  dtmb = -0.1d0*sep/dsepm
                  dtm = min(dtmb,dtm)
              END IF
          END IF
*
*       Adopt conservative time interval (expansion by factor 2 allowed).
          km0 = dtm*1.0d+03/tb
          zmu0 = body(j1)*body(j2)/body(i)
          hi = h(ipair)
          go to 50
      END IF
*
*       Specify period in yrs and N-body units.
   10 tk = semi*sqrt(semi/body(i))
      tb = yrs*tk
      tk = twopi*tk
*
*       Save energy, reduced mass & semi-major axis for corrections.
*
      hi = h(ipair)
      zmu0 = body(j1)*body(j2)/body(i)
      semi0 = semi
      ecc2 = ecc*ecc
*
*       Obtain Eddington limit during one orbit.
      dme = 2.08d-03*eddfac*(1.d0/(1.d0 + zpars(11)))*rad(2)*tb
      supedd = .false.
      novae = .false.
      disk = .false.
      rdmin = 0.0425d0*sep*(q(2)*(1.d0+q(2)))**(1.d0/4.d0)
      IF(rdmin.GT.rad(2)) disk = .true.
*
* Adopt Kelvin-Helmholz time from the modified classical expression.
*
      j = j1
      DO 15 k = 1,2
         tkh(k) = 1.0d+07*mass(k)/(rad(k)*zlmsty(j))
         IF(kstar(j).LE.1.OR.kstar(j).EQ.7.OR.kstar(j).GE.10)THEN
            tkh(k) = tkh(k)*mass(k)
         ELSE
            tkh(k) = tkh(k)*(mass(k) - massc(k))
         END IF
         j = j2
   15 CONTINUE
*
* Define dynamical timescale for the primary.
      tdyn = 5.05d-05*sqrt(rad(1)**3/mass(1))
* Identify special cases.
*
      IF(kstar(j1).EQ.2)THEN
         qc = 4.d0
      ELSEIF(kstar(j1).EQ.3.OR.kstar(j1).EQ.5.OR.kstar(j1).EQ.6)THEN
*        QC = (1.67D0-ZPARS(7)+2.D0*(MASSC(1)/MASS(1))**5)/2.13D0
*
* Alternatively use condition of Hjellming & Webbink, 1987, ApJ, 318, 794.
         qc = 0.362d0 + 1.d0/(3.d0*(1.d0 - massc(1)/mass(1)))
      ELSEIF(kstar(j1).EQ.8.OR.kstar(j1).EQ.9)THEN
         qc = 0.784d0
      ELSE
         qc = 3.d0
      ENDIF
*
      IF(kstar(j1).EQ.0.AND.q(1).GT.0.695)THEN
*
* This will be dynamical mass transfer of a similar nature to common-envelope
* evolution. The result is always a single star.
*
         taum = sqrt(tkh(1)*tdyn)
         dm1 = mass(1)
         IF (kstar(j2).LE.1)THEN
*
* Restrict accretion to thermal timescale of secondary.
*
            dm2 = taum/tkh(2)*dm1
            mass(2) = mass(2) + dm2
*
* Rejuvenate the main sequence star.
*
            mass0(2) = mass(2)
            CALL star(kstar(j2),mass0(2),mass(2),tm,tn,
     &                                 tscls,lums,gb,zpars)
* If the star has no convective core then the effective age decreases,
* otherwise it will become younger still.
            IF(mass(2).LT.0.35.OR.mass(2).GT.1.25)THEN
               aj(2) = tm/tms(2)*aj(2)*(mass(2) - dm2)/mass(2)
            ELSE
               aj(2) = tm/tms(2)*aj(2)
            END IF
            epoch(j2) = tphys0 - aj(2)
*
*       Check for blue straggler formation (TM < TPHYS).
            IF(tm.LT.tphys)THEN
               WRITE(6,16)name(j2),mass(2),tm,tphys,aj(2)
   16          FORMAT(' NEW BS (ROCHE):  NAM M3 TM TP AGE ',
     &                                   i6,f6.2,3f8.1)
            ENDIF
*
         ELSEIF(kstar(j2).LE.6)THEN
*
* Add all the material to the giant's envelope.
*
            dm2 = dm1
            mass(2) = mass(2) + dm2
            IF(kstar(j2).EQ.2)THEN
               mass0(2) = mass(2)
               CALL star(kstar(j2),mass0(2),mass(2),tm,tn,
     &                                    tscls,lums,gb,zpars)
               aj(2) = tm + tscls(1)*(aj(2)-tms(2))/tbgb(2)
               epoch(j2) = tphys0 - aj(2)
            ENDIF
         ELSEIF(kstar(j2).LE.12)THEN
*
* Form a new giant envelope (skip matrix test for preMS evolution).
*
            dm2 = dm1
            IF (min(kstar(j1),kstar(j2)).LT.0) THEN
               kst = 0
            ELSE
               kst = ktype(kstar(j1),kstar(j2))
            END IF
            IF(kst.GT.100) kst = kst - 100
            IF(kst.EQ.4)THEN
               aj(2) = aj(2)/tms(2)
               massc(2) = mass(2)
            ENDIF
*
* Check for planets or low-mass WDs.
*
            IF((kstar(j2).EQ.10.AND.mass(2).LT.0.05d0).OR.
     &         (kstar(j2).GE.11.AND.mass(2).LT.0.5d0))THEN
               kst = kstar(j1)
               mass(1) = mass(2) + dm2
               mass(2) = 0.d0
            ELSE
               mass(2) = mass(2) + dm2
               CALL gntage(massc(2),mass(2),kst,zpars,mass0(2),aj(2))
               epoch(j2) = tphys0 - aj(2)
            ENDIF
            kstar(j2) = kst
         ELSE
*
* The neutron star or black hole simply accretes at the Eddington rate.
*
            dm2 = min(dme*taum/tb,dm1)
            IF(dm2.LT.dm1) supedd = .true.
            mass(2) = mass(2) + dm2
         ENDIF
         coals = .true.
         IF(mass(2).GT.0.d0) mass(1) = 0.d0
         kw1 = kstar(j2)
         goto 60
      ELSEIF(((abs(abs(2*kstar(j1)-11)-3).EQ.2.OR.kstar(j1).EQ.9).
     &        and.(q(1).GT.qc.OR.radx(1).LE.radc(1))).OR.
     &        (kstar(j1).EQ.2.AND.q(1).GT.qc).OR.
     &        (kstar(j1).EQ.4.AND.q(1).GT.qc))THEN
*
* Common-envelope evolution for critical giants and subgiants with extreme q.
*
         m1ce = mass(1)
         m2ce = mass(2)
         WRITE (6,20)  mass, kstar(j1), kstar(j2), rad, rol, sep
   20    FORMAT (' NEW CE    M K R RL A ',2f7.3,2i3,4f9.3,f9.3)
         kw1 = kstar(j1)
         kw2 = kstar(j2)
         CALL comenv(mass0(1),mass(1),massc(1),aj(1),jspin(1),kw1,
     &               mass0(2),mass(2),massc(2),aj(2),jspin(2),kw2,
     &               ecc,sep,coals)
         WRITE (6,25)  mass, kw1, kw2, sep
   25    FORMAT (' END CE    M K A ',2f7.3,2i3,f9.3)
*
* Next step should be made without changing the time.
*
         dtm = 0.d0
         epoch(j1) = tphys0 - aj(1)
         IF(coals)THEN
            goto 60
         ELSE
            IF(kstar(j1).LT.13.AND.(kw1.EQ.13.OR.kw1.EQ.14))THEN
               jkick = j1
               kwk = kw1
            ENDIF
            IF (kstar(j2).LT.13.AND.(kw2.EQ.13.OR.kw2.EQ.14))THEN
               jkick = j2
               kwk = kw2
            ENDIF
            IF(kstar(j1).LT.10.AND.kw1.GE.10.AND.kw1.LE.12)THEN
               jkick = j1
               kwk = kw1
            ENDIF
            IF(kstar(j2).LT.10.AND.kw2.GE.10.AND.kw2.LE.12)THEN
               jkick = j2
               kwk = kw2
            ENDIF
            kstar(j1) = kw1
            kstar(j2) = kw2
            semi = sep/su
            dm1 = m1ce - mass(1)
            dm2 = mass(2) - m2ce
            dm22 = dm2
            dms(1) = 0.d0
            dms(2) = dm1 - dm2
*           jkick = 0
         ENDIF
         epoch(j2) = tphys0 - aj(2)
*
      ELSEIF(kstar(j1).GE.10.AND.kstar(j1).LE.12.AND.q(1).GE.0.628)THEN
*
* Dynamic transfer from a white dwarf.  Secondary will have KW > 9.
*
         taum = sqrt(tkh(1)*tdyn)
         dm1 = mass(1)
         IF(eddfac.LT.10.d0)THEN
            dm2 = min(dme*taum/tb,dm1)
            IF(dm2.LT.dm1) supedd = .true.
         ELSE
            dm2 = dm1
         ENDIF
         mass(2) = mass(2) + dm2
         IF(kstar(j1).EQ.10.AND.kstar(j2).EQ.10)THEN
*
* Assume the energy released by ignition of the triple-alpha reaction
* is enough to destroy the star.
*
            kstar(j2) = 15
            mass(2) = 0.d0
         ELSEIF(kstar(j1).EQ.10.OR.kstar(j2).EQ.10)THEN
*
* Should be helium overflowing onto a CO or ONe core in which case the
* helium swells up to form a giant envelope so a HeGB star is formed.
* Allowance for the rare case of CO or ONe flowing onto He is made.
*
            kstar(j2) = 9
            IF(kstar(j2).EQ.10) massc(j2) = dm2
            CALL gntage(massc(2),mass(2),kstar(j2),zpars,mass0(2),aj(2))
            epoch(j2) = tphys0 - aj(2)
         ELSEIF(kstar(j2).LE.12)THEN
            mass0(2) = mass(2)
*           IF(KSTAR(J1).EQ.12.AND.KSTAR(J2).EQ.11)THEN
*
* Mixture of ONe and CO
* or rapid accretion of CO on to CO 
* will result in an ONe product.
*
               kstar(j2) = 12
*           ENDIF
         ENDIF
         mass(1) = 0.d0
         kw1 = kstar(j2)
*       See whether Chandrasekhar's mass is exceeded (SN).
         IF(kstar(j2).LE.11.AND.mass(2).GT.mch)THEN
            WRITE (6,27)  name(j1), name(j2), kstar(j2), mass(2)
   27       FORMAT (' ROCHE SN    NAM K2* M2 ',2i6,i4,f6.2)
            mass(2) = 0.d0
            nas = nas + 1
            kw1 = 15
         ENDIF
         coals = .true.
         goto 60
*
      ELSEIF(kstar(j1).EQ.13)THEN
*
* Gamma ray burster?
*
         dm1 = mass(1)
         mass(1) = 0.d0
         kw1 = 14
         dm2 = dm1
         mass(2) = mass(2) + dm2
         coals = .true.
         ngb = ngb + 1
         goto 60
      ELSEIF(kstar(j1).EQ.14)THEN
*
* Both stars are black holes.  Let them merge quietly.
*
         dm1 = mass(1)
         mass(1) = 0.d0
         kw1 = kstar(j2)
         dm2 = dm1
         mass(2) = mass(2) + dm2
         coals = .true.
         goto 60
      ELSE
*
*       Form mass transfer for one Kepler orbit.
*
         dm1 = 3.0d-06*tb*(log(rad(1)/rol(1))**3)*
     &         min(mass(1),5.d0)**2
         IF(kstar(j1).EQ.2)THEN
            mew = (mass(1) - massc(1))/mass(1)
            dm1 = max(mew,0.01d0)*dm1
         ELSEIF(kstar(j1).GE.10)THEN
            dm1 = dm1*1.0d+03/max(rad(1),1.0d-04)
         ENDIF
         kst = kstar(j2)
*
* Limit mass transfer to the thermal rate for remaining giant-like stars and
* to the dynamical rate for all others.
*
         IF(kstar(j1).GE.2.AND.kstar(j1).LE.9.AND.kstar(j1).NE.7)THEN
            dm1 = min(dm1,mass(1)*tb/tkh(1))
         ELSEIF(rad(1).GT.10.d0*rol(1).OR.(kstar(j1).LE.1.AND.
     &          kstar(j2).LE.1.AND.q(1).GT.qc))THEN
*
* Although zeta_ad is formally infinite we do not expect a star to overfill
* its Roche lobe by more than a factor of 10 before the stars merge.
*
            WRITE (6,28)  sep, dm1, rad(1), rol(1)
   28       FORMAT (' OVERFILL    A DM1 RAD ROL ',f7.1,1p,3e10.2)
*
*       Form new star by combining the components inelastically.
            iqcoll = 3
            ch5 = ' OVER'
            goto 160
         ELSE
            dm1 = min(dm1,mass(1)*tb/tdyn)
         ENDIF
*
* Calculate wind mass loss from the stars during one orbit.
*
         vorb2 = acc1*(mass(1)+mass(2))/sep
         ivsqm = 1.d0/sqrt(1.d0-ecc*ecc)
         j = j1
         DO 30 k = 1,2
            rlperi = rol(k)*(1.d0-ecc)
            dmr(k) = mlwind(kstar(j),zlmsty(j),radx(k),
     &                      mass(k),massc(k),rlperi,zmet)
            vwind2 = 2.d0*beta*acc1*mass(k)/radx(k)
            omv2 = (1.d0 + vorb2/vwind2)**(3.d0/2.d0)
            dma(3-k) = ivsqm*acc2*dmr(k)*((acc1*mass(3-k)/vwind2)**2)/
     &                 (2.d0*sep*sep*omv2)
            dma(3-k) = min(dma(3-k),dmr(k))
            j = j2
 30      CONTINUE
*
         DO 32 k = 1,2
            dms(k) = (dmr(k)-dma(k))*tb
 32      CONTINUE
*
*       Increase time-scale to relative mass loss of 0.5% but < c.m. step.
         km = 5.0d-03/max(abs(dm1+dms(1))/mass(1),dms(2)/mass(2))
         km = min(km,2.d0*km0)
         km0 = km
         dtm = km*tb/1.0d+06
         dt1 = dtm/tstar
*
* Take the stellar evolution timestep into account but don't let it
* be overly restrictive for long lived phases.
*
         IF(iter.LE.100) dtm = min(dtm,dtmi(1),dtmi(2))
         dtm = min(dtm,(time-tev0(i))*tstar)
         dtm = max(dtm,1.0d-10)
         km = dtm*1.0d+06/tb
*
*       Decide between accreted mass by secondary and/or system mass loss.
*
         taum = mass(2)/dm1*tb
         IF(kstar(j2).LE.2.OR.kstar(j2).EQ.4)THEN
*
* Limit mass loss according to the thermal timescale of the secondary.
*
            dm2 = min(1.d0,10.d0*taum/tkh(2))*dm1
         ELSEIF(kstar(j2).GE.7.AND.kstar(j2).LE.9)THEN
*
* Naked helium star secondary swells up to a core helium burning star
* or SAGB star unless the primary is also a helium star.
*
            IF(kstar(j1).GE.7)THEN
               dm2 = min(1.d0,10.d0*taum/tkh(2))*dm1
            ELSE
               dm2 = dm1
               dmchk = dm2 - 1.05d0*dms(2)
               IF(dmchk.GT.0.d0.AND.dm2/mass(2).GT.1.0d-04)THEN
                  kst = min(6,2*kstar(j2)-10)
                  IF(kst.EQ.4)THEN
                     aj(2) = aj(2)/tms(2)
                     mcx = mass(2)
                  ELSE
                     mcx = massc(2)
                  ENDIF
                  m2 = mass(2) + km*(dm2 - dms(2))
                  CALL gntage(mcx,m2,kst,zpars,mass0(2),aj(2))
                  epoch(j2) = tphys0 + dtm - aj(2)
               ENDIF
            ENDIF
         ELSEIF(kstar(j1).LE.6.AND.
     &           (kstar(j2).GE.10.AND.kstar(j2).LE.12))THEN
*
* White dwarf secondary.
*
            IF(dm1/tb.LT.2.71d-07)THEN
               IF(dm1/tb.LT.1.03d-07)THEN
*
* Accrete until a nova explosion blows away most of the accreted material.
*
                  novae = .true.
                  dm2 = min(dm1,dme)
                  IF(dm2.LT.dm1) supedd = .true.
                  dm22 = epsnov*dm2
               ELSE
*
* Steady burning at the surface.
*
                  dm2 = dm1
               ENDIF
*
            ELSE
*
* Make a new giant envelope.
*
               dm2 = dm1
*
* Check for planets or low-mass WDs.
*
               IF((kstar(j2).EQ.10.AND.mass(2).LT.0.05d0).OR.
     &            (kstar(j2).GE.11.AND.mass(2).LT.0.5d0))THEN
                  kst = kstar(j2)
               ELSE
                  kst = min(6,3*kstar(j2)-27)
                  m2 = mass(2) + km*(dm2 - dms(2))
                  CALL gntage(massc(2),m2,kst,zpars,mass0(2),aj(2))
                  epoch(j2) = tphys0 + dtm - aj(2)
               ENDIF
            ENDIF
         ELSEIF(kstar(j2).GE.10)THEN
*
* Impose the Eddington limit.
*
            dm2 = min(dm1,dme)
            IF(dm2.LT.dm1) supedd = .true.
         ELSE
*
* We have a giant whose envelope can absorb any transferred material.
*
            dm2 = dm1
         ENDIF
         IF(.NOT.novae) dm22 = dm2
*
         IF(kst.GE.10.AND.kst.LE.12)THEN
            mt2 = mass(2) + km*(dm22 - dms(2))
            IF(kstar(j1).LE.10.AND.kst.EQ.10.AND.mt2.GE.0.7)THEN
*
* HeWD can only accrete helium-rich material up to a mass of 0.7 when
* it is destroyed in a possible Type 1a SN.
*
               mass(1) = mass(1) - km*(dm1 + dms(1))
               mass(2) = 0.d0
               kw1 = kstar(j1)
               coals = .true.
               goto 60
            ELSEIF(kstar(j1).LE.10.AND.kst.GE.11)THEN
*
* CO and ONeWDs accrete helium-rich material until the accumulated
* material exceeds a mass of 0.15 when it ignites. For a COWD with
* mass less than 0.95 the system will be destroyed in a possible
* Type 1a SN. COWDs with mass greater than 0.95 and ONeWDs will survive
* with all the material converted to ONe.
*
** Now changed to an ELD for all COWDs when 0.15 accreted (JH 11/01/00).
*
               IF((mt2-mass0(2)).GE.0.15)THEN
                  IF(kst.EQ.11)THEN
                     mass(1) = mass(1) - km*(dm1 + dms(1))
                     mass(2) = 0.d0
                     kw1 = kstar(j1)
                     coals = .true.
                     goto 60
                  ENDIF
                  mass0(2) = mt2
               ENDIF
            ELSE
               mass0(2) = mt2
            ENDIF
*
* If the Chandrasekhar limit is exceeded for a white dwarf then destroy
* the white dwarf in a supernova. If the WD is ONe then a neutron star
* will survive the supernova and we let HRDIAG take care of this when
* the stars are next updated.
* However, first check if CO WD is accreting at a fast enough rate to 
* convert to ONe (C. Tout 21/11/08). 
*
            IF(kst.EQ.11.AND.dm22.GT.0.4d0*dme/eddfac) kst = 12
            IF((kst.EQ.10.OR.kst.EQ.11).AND.mt2.GE.mch)THEN
               dm1 = mch - mass(2) + km*dms(2)
               mass(1) = mass(1) - dm1 - km*dms(1)
               mass(2) = 0.d0
               kw1 = kstar(j1)
               coals = .true.
               goto 60
            ENDIF
         ENDIF
*
*       Modify time-step & mass loss terms by speed-up factor.
         dm1 = km*dm1
         dm2 = km*dm2
         dm22 = km*dm22
         dme = km*dme
*       Save last Roche step (yrs).
         dtm0 = dtm*1.0d+06
*
* Transform to angular momentum for calculation of orbital changes.
*
         oorb = twopi/tb
         jorb = mass(1)*mass(2)/(mass(1)+mass(2))
     &          *sqrt(1.d0-ecc*ecc)*sep*sep*oorb
*
* Calculate orbital angular momentum change due to system mass loss.
*
         djorb = ((dmr(1)+q(1)*dma(1))*mass(2)*mass(2) +
     &            (dmr(2)+q(2)*dma(2))*mass(1)*mass(1))*
     &           sep*sep*oorb/(mass(1)+mass(2))**2
         djorb = djorb*dtm0
*
* For super-Eddington mass transfer rates, for gamm1 = -2.0,
* and for novae systems, assume that material is lost from
* the system as if a wind from the secondary.
* If gamm1 = -1.0 then assume the lost material carries with it
* the specific angular momentum of the primary and for all
* gamm1 > 0.0 assume that it takes away a fraction gamma of
* the orbital angular momentum.
*
*
         IF(supedd.OR.novae.OR.gamm1.LT.-1.5d0)THEN
            djorb = djorb + (dm1 - dm22)*mass(1)*mass(1)*
     &              sep*sep*oorb/(mass(1)+mass(2))**2
         ELSEIF(gamm1.GE.0.d0)THEN
            djorb = djorb + gamm1*(dm1 - dm2)*sep*sep*oorb
         ELSE
            djorb = djorb + (dm1 - dm2)*mass(2)*mass(2)*
     &              sep*sep*oorb/(mass(1)+mass(2))**2
         ENDIF
*
* For very close systems include gravitational radiation.
*
         IF(sep.LE.10.0)THEN
            CALL grrad(mass(1),mass(2),sep,ecc,jorb,djgr,delet1)
            djgr = djgr*dtm0
            djorb = djorb + djgr
            igr = igr + 1
*           IF(IGR.LT.10.OR.ABS(DJGR)/JORB.GT.0.001)THEN
            IF(sep.LT.1.05*(rad(1) + rad(2)))THEN
            IF (mod(igr,10).EQ.0) WRITE (6,45)  mass, sep, djgr, dtm
   45       FORMAT (' GR BRAKE    M1 M2 SEP DJ DTM ',2f7.3,1p,3e10.2)
            ENDIF
         ENDIF
*
         dms(1) = km*dms(1)
         IF(kstar(j1).LT.10) dms(1) = min(dms(1),mass(1) - massc(1))
         dms(2) = km*dms(2)
         IF(kstar(j2).LT.10) dms(2) = min(dms(2),mass(2) - massc(2))
*
* Calculate spin changes owing to stellar winds and magnetic braking.
*
         djspin(1) = (2.d0/3.d0)*(dmr(1)*radx(1)*radx(1)*ospin(1) -
     &                xi*dma(1)*radx(2)*radx(2)*ospin(2))*dtm0
         djspin(2) = (2.d0/3.d0)*(dmr(2)*radx(2)*radx(2)*ospin(2) -
     &                xi*dma(2)*radx(1)*radx(1)*ospin(1))*dtm0
         IF(kstar(j1).LT.10.AND.mass(1).GE.0.35)THEN
            CALL magbrk(kstar(j1),mass(1),menv(1),rad(1),ospin(1),djmb)
            djspin(1) = djspin(1) + djmb*dtm0
            imb = imb + 1
            IF(jspin(1).GT.0.d0)THEN
               IF(imb.LT.2.OR.abs(djmb)/jspin(1).GT.0.03)THEN
                  WRITE (6,40)  mass, sep, djmb, dtm
   40             FORMAT (' MB BRAKE    M1 M2 SEP DJSPN DTM ',
     &                                  2f7.3,1p,3e10.2)
               ENDIF
            ENDIF
         ENDIF
         IF(kstar(j2).LT.10.AND.mass(2).GE.0.35)THEN
            CALL magbrk(kstar(j2),mass(2),menv(2),rad(2),ospin(2),djmb)
            djspin(2) = djspin(2) + djmb*dtm0
            imb = imb + 1
            IF(jspin(2).GT.0.d0)THEN
               IF(imb.LT.10.OR.abs(djmb)/jspin(2).GT.0.001)THEN
                  WRITE (6,40)  mass, sep, djmb, dtm
               ENDIF
            ENDIF
         ENDIF
*
* Adjust the spin angular momentum of each star owing to mass transfer
* and conserve total angular momentum.
*
         djt = dm1*radx(1)*radx(1)*ospin(1)
         djspin(1) = djspin(1) + djt
         djorb = djorb - djt
         IF(disk)THEN
*
* Alter spin of the degenerate secondary by assuming that material
* falls onto the star from the inner edge of a Keplerian accretion
* disk and that the system is in a steady state.
*
            djt = dm2*twopi*aursun*sqrt(aursun*mass(2)*radx(2))
            djspin(2) = djspin(2) - djt
            djorb = djorb + djt
*
         ELSE
*
* No accretion disk.
* Calculate the angular momentum of the transferred material by
* using the radius of the disk (see Ulrich & Burger) that would
* have formed if allowed.
*
            rdisk = 1.7d0*rdmin
            djt = dm2*twopi*aursun*sqrt(aursun*mass(2)*rdisk)
            djspin(2) = djspin(2) - djt
            djorb = djorb + djt
*
         ENDIF
*
* Adjust the secondary spin if a nova eruption has occurred.
*
         IF(novae)THEN
            djt = (dm2 - dm22)*radx(2)*radx(2)*ospin(2)
            djspin(2) = djspin(2) + djt
         ENDIF
*
* Calculate any changes owing to tidal evolution (includes eccentricity).
*
         j = j1
         delet = 0.d0
         DO 500 k = 1,2
            IF((kstar(j).LE.9.AND.rad(k).GE.0.01d0*rol(k)).OR.
     &         (kstar(j).GE.10.AND.j.EQ.j1))THEN
*
               CALL bsetid(kstar(j),mass(k),massc(k),menv(k),radx(k),
     &                     radc(k),renv(k),zlmsty(j),ospin(k),k2str(k),
     &                     q(3-k),sep,ecc,oorb,delet1,dspin,eqspin,djt)
*
               delet = delet + delet1*dtm0
               IF(dtm0.GT.0.d0.AND.abs(dspin).GE.tiny)THEN
                  dspin0 = dspin
                  IF(dspin.GE.0.d0)THEN
                     dspin = min(dtm0*dspin,eqspin-ospin(k))/dtm0
                  ELSE
                     dspin = max(dtm0*dspin,eqspin-ospin(k))/dtm0
                  ENDIF
                  djt = djt*dspin/dspin0
               ELSE
                  djt = 0.d0
               ENDIF
               djorb = djorb + djt*dtm0
               djspin(k) = djspin(k) - djt*dtm0
*
            ENDIF
*
            j = j2
 500     CONTINUE
*
* Update the initial and current masses.
*
         kstar(j2) = kst
         mass(1) = mass(1) - dm1 - dms(1)
         IF(kstar(j1).LE.1.OR.kstar(j1).EQ.7) mass0(1) = mass(1)
         mass(2) = mass(2) + dm22 - dms(2)
         IF(kstar(j2).LE.1.OR.kstar(j2).EQ.7) mass0(2) = mass(2)
*
* Update spins and ensure that the star does not spin up beyond break-up.
* For a HG star check if the initial mass can be reduced.
*
         j = j1
         DO 502 k = 1,2
            jspin(k) = max(jspin(k) - djspin(k),1.0d-10)
            ospbru = twopi*sqrt(mass(k)*aursun**3/radx(k)**3)
            jspbru = (k2str(k)*(mass(k)-massc(k))*radx(k)*radx(k) +
     &                k3*massc(k)*radc(k)*radc(k))*ospbru
            IF(jspin(k).GT.jspbru)THEN
               djorb = djorb - (jspin(k) - jspbru)
               jspin(k) = jspbru
            ENDIF
            IF(kstar(j).EQ.2.AND.mass0(k).LE.zpars(3))THEN
               m01 = mass0(k)
               mass0(k) = mass(k)
               CALL star(kstar(j),mass0(k),mass(k),tmsnew,tn,tscls,
     &                   lums,gb,zpars)
               IF(gb(9).LT.massc(k)) mass0(k) = m01
            ENDIF
            j = j2
 502     CONTINUE
*
* Combine the mass lost during transfer with the secondary's wind.
*
         dms(2) = dms(2) + dm1 - dm22
*
* Adjust the orbital angular momentum and reset the orbital separation.
*
         ecc = max(ecc - delet,0.001d0)
         jorb = max(jorb - djorb,0.d0)
         sep = (mass(1) + mass(2))*jorb*jorb/
     &         ((mass(1)*mass(2)*twopi)**2*aursun**3*(1.d0-ecc*ecc))
         sep = max(sep,rad(2))
         semi = sep/su
*
*       Produce diagnostics for different types of CV and X-ray binaries.
*        WRITE (19,38)  NAME(J1), NAME(J2), KSTAR(J1), KSTAR(J2),
*    &                   TPHYS, MASS, SEP, RAD(1)/ROL(1),
*    &                   DM1/DTM/1.0E06, DM2/DTM/1.0E06, DTM
*  38    FORMAT (' ROCHE    NAM K* T M A R/RL MD12 DT ',
*    &                       2I6,2I4,F9.2,2F6.1,F7.2,F6.2,1P,3E9.1)
*        CALL flush(19)
*
*       Update KS variables to SEMI1 while maintaining low eccentricity.
*        SEMI1 = SEMI
*        CALL decay(IPAIR,SEMI1)
*
      ENDIF
*
* Rejuvenate secondary and age primary if they are on the main sequence.
*
      IF(kstar(j1).LE.2.OR.kstar(j1).EQ.7)THEN
         CALL star(kstar(j1),mass0(1),mass(1),tm,tn,
     &             tscls,lums,gb,zpars)
         IF(kstar(j1).EQ.2)THEN
            aj(1) = tm + (tscls(1)-tm)*(aj(1)-tms(1))/(tbgb(1)-tms(1))
         ELSE
            aj(1) = tm/tms(1)*aj(1)
         ENDIF
         epoch(j1) = tphys0 - aj(1)
      ENDIF
*
      IF(kstar(j2).LE.2.OR.kstar(j2).EQ.7)THEN
         CALL star(kstar(j2),mass0(2),mass(2),tm,tn,
     &             tscls,lums,gb,zpars)
         IF(kstar(j2).EQ.2)THEN
            aj(2) = tm + (tscls(1)-tm)*(aj(2)-tms(2))/(tbgb(2)-tms(2))
         ELSEIF((mass(2).LT.0.35.OR.mass(2).GT.1.25).AND.
     &      kstar(j2).NE.7)THEN
            aj(2) = tm/tms(2)*aj(2)*(mass(2) - dm22)/mass(2)
         ELSE
            aj(2) = tm/tms(2)*aj(2)
         END IF
         epoch(j2) = tphys0 - aj(2)
*
*       Check for blue straggler formation (TM < TPHYS).
         IF(kstar(j2).LE.1)THEN
            IF(tms(2).GT.tphys.AND.tm.LT.tphys)THEN
               nbr = nbr + 1
               WRITE(6,48)name(j2),mass(2),dm2,semi*su,tm,tphys,
     &                    kstar(j1)
   48          FORMAT (' NEW BR    NAM M2 DM2 A TM TP KW1 ',
     &                             i6,2f7.3,3f8.1,i4)
            ENDIF
         ENDIF
*
      ENDIF
      tphys0 = tphys0 + dtm
      tphys = tphys + dtm
*
*       Accumulate the net system mass loss.
   50 zmsy = zmsy + dms(1) + dms(2)
*
*       Update masses.
      body(j1) = mass(1)/zmbar
      body(j2) = mass(2)/zmbar
      body0(j1) = mass0(1)/zmbar
      body0(j2) = mass0(2)/zmbar
*
*       Update mass loss & binary mass and re-define binding energy.
      dmt = dms(1) + dms(2)
      dmt = dmt/zmbar
      body(i) = body(i) - dmt
      h(ipair) = -0.5d0*body(i)/semi
*
*       Include mass loss from system.
      dt2 = 1.d0
      IF (dmt.GT.0.0) THEN
          dt2 = min(dt2,0.02d0*body(j1)*dtm/(dmt*tstar))
*       Modify KS variables at constant eccentricity (expand or contract).
*         CALL expand(IPAIR,SEMI1)
          CALL expand(ipair,semi0)
*       Resolve KS components.
          CALL resolv(ipair,1)
*       Copy neighbour list for force corrections (no need!).
*         IPHASE = 0
*         NNB = LIST(1,I)
*         DO 55 L = 1,NNB+1
*             ILIST(L) = LIST(L,I)
*  55     CONTINUE
*       Correct neighbour forces and obtain the corresponding energy loss.
          IF (dmt*smu.LT.0.05) THEN
              CALL ficorr(i,dmt)
          ELSE
              CALL fcorr(i,dmt,0)
          END IF
          zmass = zmass - dmt
      ELSE IF (abs(semi - semi0).GT.1.0d-10*semi) THEN
*       Modify KS variables at constant eccentricity to treat change in SEMI.
          CALL expand(ipair,semi0)
      END IF
*
*       Subtract energy correction terms due to change in relative motion.
      zmu = body(j1)*body(j2)/body(i)
      emdot = emdot + zmu0*hi - zmu*h(ipair)
      egrav = egrav + zmu0*hi - zmu*h(ipair)
      de = zmu0*hi - zmu*h(ipair)
*
*       Increase event counter and update mass-loss.
      nroche = nroche + 1
      zmro = zmro + dm1
*
      tb = yrs*semi*sqrt(semi/body(i))
      tk = days*tb/yrs
*
*       Check termination by coalescence.
   60 IF(coals)THEN
          tphys = (time+toff)*tstar
          body0(j1) = mass0(1)/zmbar
          body0(j2) = mass0(2)/zmbar
          spin(j1) = jspin(1)/spnfac
          spin(j2) = jspin(2)/spnfac
          kw2 = 15
*       Ensure zero BODY0 for secondary to avoid confusion in COAL
*       unless both stars heve zero MASS.
          ii = j2
          k = 2
          IF(mass(1).EQ.0.0.AND.mass(2).NE.0.0)THEN
             ii = j1
             k = 1
          ENDIF
          body0(ii) = 0.d0
          tev(j1) = tphys0/tstar
          ch5 = ' COAL'
          tk = 0.0
          WRITE(85,95)name(j1),name(j2),kw1,kw2,kstar(i),
     &                tphys,aj(k),tk,mass0(k),mass0(3-k),
     &                mass(k),tk,zmet,ecc,tk,jspin(k),tk,ch5
*
          CALL coal(ipair,kw1,kw2,mass)
          go to 200
      ENDIF
*
* Obtain the stellar parameters for the next step (first J1, then J2).
*
      j = j1
      iq = 0
      DO 70 k = 1,2
         age = tphys0 - epoch(j)
         m01 = mass0(k)
         m1 = mass(k)
         mc = massc(k)
         kw = kstar(j)
         CALL star(kw,m01,m1,tm,tn,tscls,lums,gb,zpars)
         CALL hrdiag(m01,age,m1,tm,tn,tscls,lums,gb,zpars,
     &               r1,lum,kw,mc,rcc,m1ce,rce,k2)
         CALL trdot2(kw,age,tm,tn,tscls,dtmi(k),dtr)
         IF(k.EQ.1) dtmax = dtr/tstar
*
*       Update evolution times.
*
         tev(j) = tphys0/tstar
         tev0(j) = tev(j)
*        write(6,*)' roche star ',j,name(j),kstar(j),kw,log10(lum)
*
*       Quit Roche on any transition that involves a supernova as the mass
*       will change suddenly.
*
         IF((kstar(j).LT.13.AND.kw.GE.13).OR.
     &      (kstar(j).LT.10.AND.kw.GE.10.AND.kw.LE.12))THEN
            iq = j
            body0(j) = mass0(k)/zmbar
            body(j) = mass(k)/zmbar
            spin(j) = jspin(k)/spnfac
            kwk = kw
            jk = j
            go to 68
         ENDIF
*
* Save relevant solar quantities.
*
         aj(k) = age
         IF(kw.NE.kstar(j))THEN
            epoch(j) = tphys0 - age
            kstar(j) = kw
            mass0(k) = m01
            body0(j) = mass0(k)/zmbar
            mass(k) = m1
            body(j) = mass(k)/zmbar
         ENDIF
         radius(j) = r1/su
         rad(k) = r1
         radc(k) = rcc
         massc(k) = mc
         renv(k) = rce
         menv(k) = m1ce
         k2str(k) = k2
         zlmsty(j) = lum
         tms(k) = tm
         tbgb(k) = tscls(1)
         q(k) = mass(k)/mass(3-k)
         rol(k) = rl(q(k))*sep
         radx(k) = r1
         IF(j.EQ.j1)THEN
            radx(k) = min(radx(k),rol(k))
            radx(k) = max(rcc,radx(k))
         ENDIF
         ospin(k) = jspin(k)/(k2*radx(k)*radx(k)*(m1-mc) +
     &                        k3*rcc*rcc*mc)
         spin(j) = jspin(k)/spnfac
   68    j = j2
   70 CONTINUE
*
*       Check Roche termination due to sudden mass loss.
      tev0(i) = tev0(j1)
      IF(jkick.GT.0)THEN
*       Implement velocity kick without mass loss for NS & BH.
         kstar(i) = kstar(i) + 1
         ch5 = ' KICK'
         WRITE(85,95)name(j1),name(j2),kstar(j1),kstar(j2),kstar(i),
     &               tphys,aj(1),aj(2),mass0(1),mass0(2),
     &               mass(1),mass(2),zmet,ecc,tk,jspin(1),jspin(2),ch5
         CALL kick2(jkick)
         goto 160
      ELSE IF(iq.GT.0)THEN
         tev(i) = 1.0d+10
         tev(j1) = tphys0/tstar
         tev(j2) = tev(j1)
         kstar(i) = kstar(i) + 1
         kstar(jk) = kwk
         ch5 = ' KICK'
         WRITE(85,95)name(j1),name(j2),kstar(j1),kstar(j2),kstar(i),
     &               tphys,aj(1),aj(2),mass0(1),mass0(2),
     &               mass(1),mass(2),zmet,ecc,tk,jspin(1),jspin(2),ch5
         WRITE (6,71)  name(iq), name(j1), name(j2), kstar(j1),
     &                 kstar(j2), kw1, mass(1), mass(2), rad(1), rol(1)
   71    FORMAT (' ROCHE TERM    NAM K* KW M R1 RL1 ',
     &                           3i6,3i4,2f7.3,2f8.2)
         goto 140
      ENDIF
*
***
      IF(iter.EQ.0.AND.inew.EQ.1)THEN
*
*       Circularize the orbit if necessary.
*       Force co-rotation of primary with orbit if necessary.
*
         tc = 1.0d+04
         tsyn = 1.0d+20
         iterb = 0
         ecc0 = ecc
         fac0 = ospin(1)/oorb
         semi0 = semi
*
  80     CONTINUE
*
         fac = ospin(1)/oorb
         IF(fac.GT.1.d0) fac = 1.d0/fac
*
         IF((ecc-0.001d0).GT.tiny.OR.fac.LT.0.98d0)THEN
* Calculate changes owing to tidal evolution.
            delet = 0.d0
            j = j1
            DO k = 1,2
               dspink(k) = 0.d0
               djtk(k) = 0.d0
               IF((kstar(j).LE.9.AND.rad(k).GE.0.01d0*rol(k)).OR.
     &            (kstar(j).GE.10.AND.j.EQ.j1))THEN
*
                  CALL bsetid(kstar(j),mass(k),massc(k),menv(k),radx(k),
     &                      radc(k),renv(k),zlmsty(j),ospin(k),k2str(k),
     &                      q(3-k),sep,ecc,oorb,delet1,dspin,eqspin,djt)
*
                  delet = delet + delet1
                  dspink(k) = dspin
                  djtk(k) = djt
               ENDIF
               j = j2
            ENDDO
* Include gravitational radiation.
            CALL grrad(mass(1),mass(2),sep,ecc,jorb,djgr,delet1)
            delet = delet + delet1
* Include magnetic braking.
            CALL magbrk(kstar(j1),mass(1),menv(1),rad(1),ospin(1),djmb)
            djspin(1) = djtk(1) + djmb
*
            ttid = 1.0d+20
            dtm0 = 1.0d+20
            dtmmin = 0.d0
            IF(ecc.GT.0.001d0)THEN
               IF(abs(delet).GT.1.0d-14)THEN
* Allow eccentricity to change by 10% at most with a minimum of 0.01.
                  dtm0 = max(0.1d0*ecc,0.0005d0)/abs(delet)
                  dtmmin = min(0.01d0,ecc-0.001d0)/abs(delet)
                  tc = ecc/abs(delet)
                  tc = tc/1.0d+06
                  ttid = tc
               ELSE
* Perform instant circularization conserving orbital angular momentum. 
                  ecc = 0.001d0
                  sep = jorb*(mass(1)+mass(2))/(mass(1)*mass(2)*oorb)
                  sep = sqrt(sep)
                  tb = (sep/aursun)*sqrt(sep/(aursun*(mass(1)+mass(2))))
                  tk = days*tb/yrs
                  oorb = twopi/tb
               ENDIF
            ELSE
               delet = 0.d0
            ENDIF
            IF(abs(djspin(1)).GT.1.0d-14)THEN
               IF(abs(dspink(1)).GT.1.0d-20)THEN
                  tsyn = abs(oorb-ospin(1))/abs(dspink(1))
                  dtm0 = min(dtm0,0.3d0*tsyn)
                  tsyn = tsyn/1.0d+06
                  ttid = min(ttid,tsyn)
               ELSE
* Allow jspin to change by 50% at most.
                  dtm0 = min(dtm0,0.5d0*jspin(1)/abs(djspin(1)))
               ENDIF
            ENDIF
*
            IF(iterb.EQ.0.AND.nwarn.LT.50)THEN
               WRITE(6,710)name(j1),kstar(j1),kstar(j2),ecc0,fac0,
     &                     tc,tsyn
  710          FORMAT(' CIRC & SYNCH NM K* ECC0 SPIN1/OORB TC TSYN ',
     &                                i7,2i4,f7.4,f9.3,1p,2e10.2)
               nwarn = nwarn + 1
            ENDIF
*
            IF(ttid.LT.min(tphys,10.d0))THEN
               itro = 0
               dtm0 = max(dtm0,dtmmin)
               delet = delet*dtm0
   81          djorb = djgr*dtm0
* Include magnetic braking for both stars.
               djspin(1) = djmb*dtm0
               CALL magbrk(kstar(j2),mass(2),menv(2),rad(2),ospin(2),
     &                     djmb)
               djspin(2) = djmb*dtm0
* Check equilibrium spin and add tidal changes.
               DO k = 1,2
                  IF(dspink(k).GT.1.0d-10)THEN
                     dspin0 = min(dtm0*dspink(k),eqspin-ospin(k))/dtm0
                     djtk(k) = dspin0*djtk(k)/dspink(k)
                  ELSEIF(dspink(k).LT.-1.0d-10)THEN
                     dspin0 = max(dtm0*dspink(k),eqspin-ospin(k))/dtm0
                     djtk(k) = dspin0*djtk(k)/dspink(k)
                  ELSE
                     dspink(k) = 0.d0
                  ENDIF
                  djorb = djorb + djtk(k)*dtm0
                  djspin(k) = djspin(k) - djtk(k)*dtm0
               ENDDO
* Update spins and ensure that the star does not spin up beyond break-up.
               DO k = 1,2
                  jspin(k) = max(jspin(k) - djspin(k),1.0d-10)
                  ospbru = twopi*sqrt(mass(k)*aursun**3/radx(k)**3)
                  jspbru = (k2str(k)*(mass(k)-massc(k))*radx(k)*radx(k)
     &                      + k3*massc(k)*radc(k)*radc(k))*ospbru
                  IF(jspin(k).GT.jspbru)THEN
                     djorb = djorb - (jspin(k) - jspbru)
                     jspin(k) = jspbru
                  ENDIF
                  ospin(k) = jspin(k)/
     &                      (k2str(k)*radx(k)*radx(k)*(mass(k)-massc(k))
     &                       + k3*radc(k)*radc(k)*massc(k))
               ENDDO
* Impose 2.5 % limit on change of JORB (Chris Tout 1/2013).
               IF (abs(djorb).GT.0.1*jorb.AND.itro.LE.5) THEN
                   dtm0 = 0.1*jorb*dtm0/abs(djorb)
                   delet = 0.1*jorb*delet/abs(djorb)
                   itro = itro + 1
                   IF (jorb.GT.0.0d0) WRITE (6,83)  djorb/jorb, jorb
   83              FORMAT (' ROCHE LIMIT    DJO/JORB JORB ',1p,2e10.2)
                   go to 81
               ENDIF
* Adjust the eccentricity and orbital angular momentum.
               ecc = max(ecc - delet,0.001d0)
               jorb = max(jorb - djorb,0.d0)
* Reset the orbital separation and Roche-lobe radii.
               sep = (mass(1) + mass(2))*jorb*jorb/
     &             ((mass(1)*mass(2)*twopi)**2*aursun**3*(1.d0-ecc*ecc))
               sep = max(sep,rad(2))
* Perform coalescence check after shrinkage.
               IF (rad(1)+rad(2).LT.sep)THEN
                  coals = .true.
                  go to 60
               ENDIF
               DO k = 1,2
                  rol(k) = rl(q(k))*sep
               ENDDO
               radx(1) = min(rad(1),rol(1))
               radx(1) = max(radc(1),radx(1))
*
               tb = (sep/aursun)*sqrt(sep/(aursun*(mass(1)+mass(2))))
               tk = days*tb/yrs
               oorb = twopi/tb
               iterb = iterb + 1
               IF(iterb.LE.100) goto 80
            ELSE
               IF(ecc.GT.0.02d0)THEN
                  WRITE(6,*)' WARNING: will not circularize'
               ENDIF
               IF(fac.LT.0.98d0.AND.kstar(j1).LT.10)THEN
                  WRITE(6,*)' WARNING: will not synchronize'
               ENDIF
               iterb = iterb + 1
            ENDIF
         ENDIF
         IF(iterb.GT.0)THEN
            WRITE(6,711)ecc,ospin(1)/oorb,sep,tk,rad(1),rol(1)
  711       FORMAT(' NEW PARAMS    ECC SPIN1/OORB A P R1 RL1 ',
     &                             f7.4,f9.3,f8.2,1p,3e10.2)
            spin(j1) = jspin(1)/spnfac
            spin(j2) = jspin(2)/spnfac
*       Set new orbital elements (ECC, SEMI & JORB have changed).
            hi = h(ipair)
            semi = sep/su
*       Specify consistent orbital energy from new semi-major axis.
            h(ipair) = -0.5d0*body(i)/semi
*       Change the KS variables at constant eccentricity.
            CALL expand(ipair,semi0)
*       Transform to pericentre (note TIME may change).
            CALL ksperi(ipair)
*       Modify the eccentricity at constant energy.
            CALL deform(ipair,ecc0,ecc)
            zmu = body(j1)*body(j2)/body(i)
            emdot = emdot + zmu*(hi - h(ipair))
            egrav = egrav + zmu*(hi - h(ipair))
         ENDIF
*
      ENDIF
*
*       Check if components should be saved for output, cf. HRPLOT
      j = j1
      IF(tev(j).GE.tplot.AND..NOT.isave)THEN
         isave = .true.
*       Check if components are already saved (temporary)
         ii = 0
         DO 700 k = 1,nrsave
            IF(name(j1).EQ.namer(1,k)) ii = name(j1)
            IF(name(j1).EQ.namer(2,k)) ii = name(j1)
  700    CONTINUE
         IF(ii.EQ.0)THEN
            nrsave = nrsave + 1
            IF(nrsave.GT.mmax)THEN
               iwarn = iwarn + 1
               IF (iwarn.LT.20) THEN
                   WRITE(6,*)' ROCHE WARNING: NRSAVE > MMAX'
               END IF
               nrsave = mmax
            ENDIF
            namer(1,nrsave) = name(j1)
            namer(2,nrsave) = name(j2)
            kstarr(1,nrsave) = kstar(i)
            kstarr(2,nrsave) = kstar(j1)
            kstarr(3,nrsave) = kstar(j2)
            cmrch(1,nrsave) = semi
            DO 701 k = 1,3
               cmrch(k+1,nrsave) = x(k,i)
               cmrch(k+4,nrsave) = xdot(k,i)
  701       CONTINUE
            cmrch(8,nrsave) = body0(j1)
            cmrch(9,nrsave) = body(j1)
            cmrch(10,nrsave) = epoch(j1)
            cmrch(11,nrsave) = body0(j2)
            cmrch(12,nrsave) = body(j2)
            cmrch(13,nrsave) = epoch(j2)
         ENDIF
      ENDIF
*
* See whether the primary still fills its Roche lobe.
*
      IF(rad(1).GT.rol(1))THEN
*
* Test for a contact system.
*
         IF (rad(2).GT.rol(2).AND.iter.GT.0) THEN
            ncont = ncont + 1
            WRITE (6,72)  mass(1), mass(2), kstar(j1), kstar(j2), sep
   72       FORMAT (' CONTACT    M1 M2 KW1 KW2 A ',2f7.3,2i3,f8.3)
*
*       Form new star by combining the components inelastically.
            iqcoll = 3
            ch5 = ' CONT'
            go to 160
         END IF
*
         IF (iwarn.LT.100.AND.rad(1).GT.2.0*rol(1)) THEN
            iwarn = iwarn + 1
            WRITE (6,73)  kstar(j1), kstar(j2), mass(1), mass(2),
     &                    rad(1), rol(1), dm1, dm2, dtm
   73       FORMAT (' BIG ROCHE    K12 M12 R1 RL1 DM12 DT ',
     &                             2i3,2f7.3,2f8.3,1p,3e10.2)
          END IF
*
*       Produce diagnostics for cases involving degenerate objects.
          IF (iter.EQ.1.AND.inew.GT.0.AND.kstar(j2).GE.10) THEN
              IF (abs(dtm).GT.1.0d-20) THEN
                  pd = days*semi*sqrt(semi/body(i))
                  WRITE (22,58)  name(j1), name(j2),
     &                           kstar(j1), kstar(j2),
     &                           mass(1), mass(2), tphys, semi*su, pd,
     &                           dm1/dtm/1.0e+06, dm2/dtm/1.0e+06
   58             FORMAT (' DEGEN ROCHE    NAM K* M TP A P M1D M2D ',
     &                             2i6,2i4,2f6.2,f8.1,2f7.1,1p,2e9.1)
                  CALL flush(22)
              END IF
          END IF
*
          iter = iter + 1
          IF(tev0(i).LT.time) goto 10
          IF (iter.EQ.1.AND.gamma(ipair).LT.gmin) go to 10
      ELSE
          WRITE (6,76)  kstar(j1), kstar(j2), mass(1),
     &                  mass(2),rad(1), rol(1), dm1, dm2, dtm
   76     FORMAT(' END ROCHE    K12 M12 R1 RL1 DM DT ',
     &                          2i3,2f7.3,2f8.2,2f8.3,1p,e10.2)
*       Check optional diagnostics for degenerate objects.
          IF(max(kstar(j1),kstar(j2)).GE.10)THEN
             IF(kz(8).GT.3)THEN
                CALL degen(ipair,ipair,4)
             ELSE
                ndd = ndd + 1
             ENDIF
          ENDIF
          IF (kstar(j1).GE.10.AND.kstar(j1).LE.12) nwd = nwd + 1
*       Define end of Roche stage.
          kstar(i) = kstar(i) + 1
          CALL trflow(ipair,dtr)
          tev(i) = time + dtr
          IF (kz(8).GT.3) THEN
              CALL binev(ipair)
          END IF
          ch5 = ' END '
          WRITE(85,95)name(j1),name(j2),kstar(j1),kstar(j2),kstar(i),
     &                tphys,aj(1),aj(2),mass0(1),mass0(2),
     &                mass(1),mass(2),zmet,ecc,tk,jspin(1),jspin(2),ch5
          CALL flush(85)
      END IF
*
*       Check standard exit condition.
      IF(iter.EQ.0) goto 150
*
  140 CONTINUE
*       Set new c.m. time only during active stage (otherwise by TRFLOW).
      IF (mod(kstar(i),2).EQ.1.AND.iq.EQ.0) THEN
          tev(i) = tphys0/tstar
*       Advance update time (extra for critical configurations).
          IF(list(1,i1).GT.0)THEN
             dtmax = min(dtmax,50.d0*dt1)
          ELSE
             dtmax = min(dtmax,10.d0*dt1)
          ENDIF
          dtmax = min(dtmax,dt2)
*         DTMAX = MIN(DTMAX,STEPX)
          tev(i) = tev(i) + max(dtmax,1.0d-07)
          tev(j1) = tev(i) + 2.0*stepx + step(i1)
          tev(j2) = tev(j1)
          tmdot = min(tmdot,tev(i))
          z = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
          ecc = sqrt(z)
          IF (ecc.GT.0.01) THEN
              qp = semi*(1.0 - ecc)
              icirc = -1
              CALL tcirc(qp,ecc,i1,i2,icirc,tc)
              WRITE (6,142)  time+toff, name(j1), kstar(j1), kstar(j2),
     &                       list(1,i1), mass(1), mass(2), ecc,
     &                       semi*su, tc
  142         FORMAT (' ECCENTRIC ROCHE    T NM K* NP M E A TC ',
     &                                     f8.1,i6,3i4,3f6.2,2f7.1)
              kstar(i) = 0
              tev(i) = 1.0d+10
          END IF
      ELSEIF(iq.EQ.0)THEN
*       Update look-up time.
          CALL trdot(j1,dtm,mass(1))
          tev(j1) = tev(j1) + dtm
          CALL trdot(j2,dtm,mass(2))
          tev(j2) = min(tev(j1),tev(j2) + dtm)
*
*       Ensure that close degenerate system gets GR check in mdot. 
          kw = min(kstar(j1),kstar(j2))
          IF(kw.GT.6) tev(j2) = min(tev(j2),time+0.1d0)
*
          tev(j1) = tev(j2)
      END IF
*
*       Restore times and step and update mass on GRAPE.
      tphys = (time+toff)*tstar
*
*       Include consistency check on mass difference.
      dm = (body(j1) + body(j2)) - body(i)
      IF(abs(dm).GT.0.0001*body(i))THEN
          WRITE (6,145)  time+toff, body(j1), body(j2), body(i)
  145     FORMAT (' DANGER!    ROCHE    T M ',f10.3,1p,3e12.4)
          body(i) = body(i) + dm
      ENDIF
*
      go to 200
*
  150 IF(iphase.GE.0) goto 140
  160 tphys = (time+toff)*tstar
*
*       Treat case of OVERFILL or CONTACT binary in the same way.
      IF(iqcoll.EQ.3)THEN
*       Predict c.m. to highest order and save KS index in COMMON.
          WRITE(85,95)name(j1),name(j2),kstar(j1),kstar(j2),kstar(i),
     &                tphys,aj(1),aj(2),mass0(1),mass0(2),
     &                mass(1),mass(2),zmet,ecc,tk,jspin(1),jspin(2),ch5
          IF(time-t0(i).LE.step(i))THEN
              CALL xvpred(i,0)
          ENDIF
          kspair = ipair
          CALL cmbody(r(ipair),2)
          iqcoll = 0
      ENDIF
*
  200 RETURN
*
      END
***