SUBROUTINE TRFLOW(IPAIR,DTR)
*
*
*       Time until Roche overflow.
*       --------------------------
*
*
      INCLUDE 'common6.h'
      REAL*8  TSCLS(20),LUMS(10),GB(10),TM,TN
      REAL*8  M0,M1,AGE,LUM,MC,RCC,Q,RL1,RL2
      REAL*8  MENV,RENV,K2
      REAL*8  AJ,RHIGH,THIGH,T,TT,RR,RX,TLOW,RLOW,DEL,DER,EPS,TOL
      PARAMETER (EPS = 1.0D-06, TOL = 1.0D-04)
      REAL*8 RL,RTMSF,RGBF,RGBDF,RAGBF,RAGBDF,RZHEF
      EXTERNAL RL,RTMSF,RGBF,RGBDF,RAGBF,RAGBDF,RZHEF
      LOGICAL ITEV
*
*
*       Set index of first component and semi-major axis.
      ITEV = .TRUE.
      IF(IPAIR.LT.0)THEN
         ITEV = .FALSE.
         IPAIR = -IPAIR
      ENDIF
      I1 = 2*IPAIR - 1
      I = N + IPAIR
      SEMI = -0.5*BODY(I)/H(IPAIR)
*
*       Determine indices for primary & secondary star (donor & accretor).
      J1 = 2*IPAIR - 1
      J2 = I1 + 1
      Q = BODY(J1)/BODY(J2)
      RL1 = RL(Q)*SEMI
*       Evaluate Roche radius for the second star.
      Q = 1.0/Q
      RL2 = RL(Q)*SEMI
*
*       Compare scaled Roche radii when choosing the primary.
      IF (RADIUS(J1)/RL1.LT.RADIUS(J2)/RL2) THEN
          J1 = I1 + 1
          J2 = I1
          RL1 = RL2
      END IF
*     IF(TEV(J1).LE.0.0)THEN
*        WRITE(6,99)J1,NAME(J1),NAME(I),KSTAR(J1),TEV(J1)+TOFF,TTOT
*99      FORMAT(' TRFLOW WARNING! ',3I6,I4,2F10.3)
*     ENDIF
*
*       Exit with large interval if semi-major axis is negative.
      IF(SEMI.LE.0.D0)THEN
          DTR = 1.0D+10
          GOTO 210
*       Assign DTR for active Roche during coasting.
      ELSEIF(KSTAR(I).GT.0.AND.MOD(KSTAR(I),2).NE.0)THEN
          IF(TEV(I).GT.9.9D+09) TEV(I) = TEV0(J1)
          DTR = TEV(I) - TIME
          GOTO 210
      ENDIF
*
*       Convert Roche radius, radius and initial & current mass to SU.
      RLS = RL1*SU
      M0 = BODY0(J1)*ZMBAR
      M1 = BODY(J1)*ZMBAR
      MC = 0.D0
      KW = KSTAR(J1)
      CALL star(KW,M0,M1,TM,TN,TSCLS,LUMS,GB,ZPARS)
      AGE = TEV0(J1)*TSTAR - EPOCH(J1)
      CALL hrdiag(M0,AGE,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &            RSJ,LUM,KW,MC,RCC,MENV,RENV,K2)
*
* If the star already fills its Roche lobe exit with DTR = 0.
*
      IF(RSJ.GE.RLS)THEN
         DTR = 0.D0
         RADIUS(J1) = RSJ/SU
         GOTO 210
      ENDIF
*
*       Exit with large interval if primary has terminated its evolution.
      IF(KSTAR(J1).GE.10)THEN
          DTR = 1.0D+10
          GOTO 210
      ENDIF
*
* If the star is on the MS then see if it will fill its Roche Lobe
* during this stage.
*
      T = 2.0D+10
      IF(KW.LE.1.OR.KW.EQ.7)THEN
         AJ = 0.99D0*TM
         CALL hrdiag(M0,AJ,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &               RHIGH,LUM,KW,MC,RCC,MENV,RENV,K2)
         IF(RHIGH.GT.RLS.AND.AJ.GT.AGE)THEN
* Find the overflow time using binary chopping.
            TLOW = AGE
            RLOW = RSJ
            THIGH = AJ
            IT = 0
   10       T = 0.5D0*(THIGH + TLOW)
            CALL hrdiag(M0,T,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &                  RR,LUM,KW,MC,RCC,MENV,RENV,K2)
            IF(RR.LT.RLS)THEN
               TLOW = T
               RLOW = RR
            ELSE
* Stop when 0.1% accuracy is achieved.
               IT = IT + 1
               IF(IT.EQ.25)THEN
                  IF(ABS(RR-RLS).GT.0.1)THEN
                     WRITE(38,*)' TRFLOW KW1: ',rr,rls,t,tm
                  ENDIF
               ENDIF
               IF(ABS(RR - RLS).LE.0.001*RR.OR.IT.GT.25) GOTO 200
               THIGH = T
               RHIGH = RR
            ENDIF
            GOTO 10
         ENDIF
      ENDIF
*
* See if the star will fill its Roche Lobe on the HG if not already
* past this stage.
*
      IF(KW.LE.2)THEN
         AJ = MIN(TN,TSCLS(1))*(1.D0-EPS)
         CALL hrdiag(M0,AJ,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &               RHIGH,LUM,KW,MC,RCC,MENV,RENV,K2)
         IF(RHIGH.GT.RLS.AND.AJ.GT.AGE)THEN
* Solve for the overflow time.
            RR = RTMSF(M0)
            T = LOG(RLS/RR)/LOG(RHIGH/RR)
            T = TM + T*(TSCLS(1) - TM)
            GOTO 200
         ENDIF
         IF(TN.LE.TSCLS(1)) GOTO 200
      ENDIF
*
* If the star is pre-helium ignition see if it will fill its Roche
* lobe before helium ignition.
*
      IF(KW.LE.3)THEN
         AJ = MIN(TN,TSCLS(2))*(1.D0-EPS)
         CALL hrdiag(M0,AJ,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &               RHIGH,LUM,KW,MC,RCC,MENV,RENV,K2)
         IF(RHIGH.GE.RLS)THEN
* Solve for the luminosity corresponding to the Roche radius and
* then for the overflow time.
            IT = 0
            LUM = LUMS(3)
 30         IT = IT + 1
            DEL = RGBF(M1,LUM) - RLS
            IF(IT.EQ.25.AND.ABS(DEL).GT.0.1)THEN
               WRITE(38,*)' TRFLOW KW3: ',rgbf(m1,lum),rls,lum,lums(3)
            ENDIF
            IF(ABS(DEL/RLS).LE.TOL.OR.IT.EQ.25) GOTO 40
            DER = RGBDF(M1,LUM)
            LUM = LUM - DEL/DER
            GOTO 30
 40         CONTINUE
            IF(LUM.LE.LUMS(6))THEN
               T = TSCLS(4) - (1.D0/((GB(5)-1.D0)*GB(1)*GB(4)))*
     &                      ((GB(4)/LUM)**((GB(5)-1.D0)/GB(5)))
            ELSE
               T = TSCLS(5) - (1.D0/((GB(6)-1.D0)*GB(1)*GB(3)))*
     &                      ((GB(3)/LUM)**((GB(6)-1.D0)/GB(6)))
            ENDIF
            GOTO 200
         ELSE
* If a low mass star has not yet filled its Roche Lobe then it will
* not do it until after the Helium Flash. Thus we don't let it go
* any further until it has actually become a type 4.
            IF(M0.LT.ZPARS(2).OR.TN.LE.TSCLS(2)) GOTO 200
         ENDIF
      ENDIF
*
* Check for overflow during the CHeB stage.
*
      IF(KW.EQ.4)THEN
*
         IF(TN.LT.(TSCLS(2)+TSCLS(3)))THEN
            T = MAX(TSCLS(2),AGE)
            AJ = T + 0.5D0*(TN - T)
            CALL hrdiag(M0,AJ,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &                  RHIGH,LUM,KW,MC,RCC,MENV,RENV,K2)
            IF(RHIGH.LT.RLS)THEN
* If the evolution is due to end during CHeB then the pertubation
* functions for small envelope mass will take effect at some point
* causing the radius to decrease. We assume that this point is after
* AJ and quit with T as an underestimate of the overflow time.
               T = AJ
               GOTO 200
            ENDIF
         ELSE
            AJ = (TSCLS(2)+TSCLS(3))*(1.D0-EPS)
            CALL hrdiag(M0,AJ,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &                  RHIGH,LUM,KW,MC,RCC,MENV,RENV,K2)
            IF(RHIGH.LT.RLS) GOTO 50
         ENDIF
* Find the overflow time using binary chopping.
         TLOW = AGE
         RLOW = RSJ
         THIGH = AJ
         IT = 0
 60      T = 0.5D0*(THIGH + TLOW)
         CALL hrdiag(M0,T,M1,TM,TN,TSCLS,LUMS,GB,ZPARS,
     &               RR,LUM,KW,MC,RCC,MENV,RENV,K2)
         IF(RR.LT.RLS)THEN
            TLOW = T
            RLOW = RR
         ELSE
* Stop when 0.1% accuracy is achieved.
            IT = IT + 1
            IF(IT.EQ.25)THEN
               IF(ABS(RR-RLS).GT.0.1)THEN
                  WRITE(38,*)' TRFLOW KW4: ',rr,rls,tscls(2),t,tscls(3)
               ENDIF
            ENDIF
            IF(ABS(RR - RLS).LE.0.001*RR.OR.IT.GT.25) GOTO 200
            THIGH = T
            RHIGH = RR
         ENDIF
         GOTO 60
 50      CONTINUE
      ENDIF
*
* The star now has only until the end of the AGB phase to fill its lobe.
*
      IF(KW.LE.6)THEN
* Solve for the luminosity corresponding to the Roche radius and
* then for the overflow time.
         IT = 0
         IF(KW.EQ.6)THEN
            LUM = LUMS(8)
         ELSE
            LUM = LUMS(7)
         ENDIF
         RR = RAGBF(M1,LUM,ZPARS(2))
         IF(RR.GT.RLS)THEN
* In this case the solution should already have been found. If it hasn't
* then most likely the envelope is small and the pertubation functions
* are taking effect.
            T = TN
            GOTO 200
         ENDIF
 70      IT = IT + 1
         RR = RAGBF(M1,LUM,ZPARS(2))
         DEL = RR - RLS
         IF(IT.EQ.25.AND.ABS(DEL).GT.0.1)THEN
            WRITE(38,*)' TRFLOW KW6: ',rr,rls,lum,lums(7)
         ENDIF
         IF(ABS(DEL/RLS).LE.TOL.OR.IT.EQ.25) GOTO 80
         DER = RAGBDF(M1,LUM,ZPARS(2))
         LUM = LUM - DEL/DER
         GOTO 70
 80      CONTINUE
         IF(LUM.LE.LUMS(8))THEN
            IF(LUM.LE.LUMS(6))THEN
               T = TSCLS(7) - (1.D0/((GB(5)-1.D0)*GB(8)*GB(4)))*
     &                      ((GB(4)/LUM)**((GB(5)-1.D0)/GB(5)))
            ELSE
               T = TSCLS(8) - (1.D0/((GB(6)-1.D0)*GB(8)*GB(3)))*
     &                      ((GB(3)/LUM)**((GB(6)-1.D0)/GB(6)))
            ENDIF
         ELSE
            IF(LUM.LE.LUMS(6))THEN
               T = TSCLS(10) - (1.D0/((GB(5)-1.D0)*GB(2)*GB(4)))*
     &                       ((GB(4)/LUM)**((GB(5)-1.D0)/GB(5)))
            ELSE
               T = TSCLS(11) - (1.D0/((GB(6)-1.D0)*GB(2)*GB(3)))*
     &                       ((GB(3)/LUM)**((GB(6)-1.D0)/GB(6)))
            ENDIF
         ENDIF
         GOTO 200
      ENDIF
*
      IF(KW.LE.9)THEN
*
* The star has until the end of the Helium GB to fill its lobe.
*
         LUM = (RLS/0.08D0)**(4.0/3.0)
         CM = 2.0D-03*M1**2.5/(2.D0 + M1**5)
         IF(CM*LUM.GE.100.D0)THEN
            DTR = 1.0D+10
            GOTO 210
         ENDIF
         RX = RZHEF(M1)
         RR = RX*(LUM/LUMS(2))**0.2 + 
     &            0.02D0*(EXP(CM*LUM) - EXP(CM*LUMS(2)))
         IF(RR.LT.RLS)THEN
* Find the overflow luminosity using binary chopping.
            TLOW = LUM
            RLOW = RR
 89         LUM = 2.0*LUM
            RR = RX*(LUM/LUMS(2))**0.2 + 
     &               0.02D0*(EXP(CM*LUM) - EXP(CM*LUMS(2)))
            IF(RR.LE.RLS) GOTO 89
            THIGH = LUM
            RHIGH = RR
            IT = 0
 90         LUM = 0.5D0*(THIGH + TLOW)
            RR = RX*(LUM/LUMS(2))**0.2 + 
     &               0.02D0*(EXP(CM*LUM) - EXP(CM*LUMS(2)))
            IF(RR.LT.RLS)THEN
               TLOW = LUM
               RLOW = RR
            ELSE
* Stop when 0.1% accuracy is achieved.
               IT = IT + 1
               IF(IT.EQ.25)THEN
                  IF(ABS(RR-RLS).GT.0.1)THEN
                     WRITE(38,*)' TRFLOW KW9: ',rr,rls,t,tm
                  ENDIF
               ENDIF
               IF(ABS(RR - RLS).LE.0.001*RR.OR.IT.GT.25) GOTO 95
               THIGH = LUM
               RHIGH = RR
            ENDIF
            GOTO 90
 95         CONTINUE
         ENDIF
         IF(LUM.LE.LUMS(6))THEN
            T = TSCLS(4) - (1.D0/((GB(5)-1.D0)*GB(8)*GB(4)))*
     &                   ((GB(4)/LUM)**((GB(5)-1.D0)/GB(5)))
         ELSE
            T = TSCLS(5) - (1.D0/((GB(6)-1.D0)*GB(8)*GB(3)))*
     &                   ((GB(3)/LUM)**((GB(6)-1.D0)/GB(6)))
         ENDIF
      ENDIF
*
 200  CONTINUE
      IF(T.GE.TN)THEN
         DTR = 1.0D+10
      ELSE
         TT = (T+EPOCH(J1))/TSTAR
         DTR = TT - TIME
         IF(DTR.LT.0.D0.AND.ITEV)THEN
            TEV(J1) = MIN(TEV(J1),TIME)
            TEV(J2) = TEV(J1)
         ENDIF
         DTR = MAX(DTR,0.D0)
         IF(TEV(J1).LT.TEV0(J1))THEN
            WRITE (6,101)J1,KSTAR(J1),TEV(J1)+TOFF,TTOT
 101        FORMAT(' TRFLOW WARNING! TEV<TEV0 ',I6,I4,2F10.3)
         ENDIF
      ENDIF
*
 210  RETURN
      END
***