      subroutine root(a,b,f,tol,flag)
      implicit none
      double precision a,b,f,tol
      integer flag
c*****************************************************************
c  subroutine root: find a root of a given function of one
c                   variable using a hybrid algorithm involving
c                   the bisection and secant methods
c
c  Author: Douglas N. Arnold, 9/96
c
c  input:
c
c  a,b      endpoints of a bracketing interval
c  f        function subroutine f(x)
c  tol      desired length of the final bracketing interval
c
c  output:
c
c  a,b      endpoints of a new bracketing interval
c  flag     0 if tolerance achieved; in this case, the bracketing
c             interval returned has size at most tol
c           1 if tolerance could not be achieved; in this case the
c             returned bracketing interval is as small as possible
c           2 if exact root is found; in this case a and b are
c             both set equal to the root
c          -1 if input data in error; in this case a and b are
c             left unchanged from their input values
c*****************************************************************
      double precision temp,c,bnew,fa,fb,fc,fbnew,delta,deltas,
     +       deltam,deltaf
      fa = f(a)
      fb = f(b)
      c = a
      fc = fa
c check for exact root in input data
      if (fa .eq. 0.d0) then
        b=a
        flag = 2
        return
      endif
      if (fb .eq. 0.d0) then
        a=b
        flag = 2
        return
      endif
c check for bad input data
      if (sign(1.d0,fa) .eq. sign(1.d0,fb) .or. tol .le. 0.d0 ) then
         flag = -1
         return
      endif
      
c main loop
   10 continue
        if (abs(fb) .gt. abs(fa)) then
          temp = b
          b = a
          a = temp
          temp = fb
          fb = fa
          fa = temp
          c = a
          fc = fa
        endif
c   check for convergence
        if (abs(a-b) .le. tol) then
          flag = 0
          return
        endif
c   compute the increment for the midpoint method based on a and b
        deltam = (a-b)/2.d0
c   if interval is as short as possible, return
        if ( b + deltam .eq. b) then
          flag = 1
          return
        endif
c   compute the increment for the secant method based on b and c
        deltaf = fc - fb
        if (deltaf .eq. 0.d0) then
          deltas = deltam
        else
          deltas = -fb*(c-b)/deltaf
        endif
c   if the secant increment seems reasonable, use that to update b
c   else, use the midpoint increment
        if (sign(1.d0,deltas) .eq. sign(1.d0,deltam)  .and. 
     +      abs(deltas) .lt. abs(deltam) ) then
c          print *, 'secant'
          delta = deltas
        else
c          print *,'midpoint'
          delta = deltam
        endif
c   if the increment is very small, move by 1/2 tol
        if (abs(delta) .lt. tol) then
c          print*,'adjusting for small increment'
          delta = 0.5*sign(1.d0,deltam)*tol
        endif
c   update
        bnew = b + delta
        fbnew = f(bnew)
        if (fbnew .eq. 0.d0) then
          a = bnew
          b = bnew
          flag = 2
          return
        endif
        c = b
        b = bnew
        fc = fb
        fb = fbnew
c    maintain a bracketing interval
        if (sign(1.d0,fb) .eq. sign(1.d0,fa)) then
          a = c
          fa = fc
        endif
      go to 10
      
      end
