Subroutine root_bs(f,a,b,eps,Root,iflag)
c====================================================================
c Program to find a single root of an equation f(x)=0
c using the Bisectional method
c written by: Alex Godunov (October 2006)
c--------------------------------------------------------------------
c input ...
c f	- function which evaluates f(x) for any x in the interval a,b
c a	- left endpoint of initial interval
c b	- right endpoint of initial interval
c eps	- desired uncertainity of the root
c
c output ...
c Root  - root of equation f(x)=0
c iflag - indicator of success
c         1 - a single root found
c         0 - no solutions for the Bisectional method
c
c Limitations: a function f(x) must change sign between a and b
c
c Max number of allowed iterations is 1000 (accuracy (b-a)/2**1000)
c====================================================================
	Real*8 f,a,b,eps,Root
	Real*8 xl,x0,xr
	Integer*4 i, iter, iflag
	Parameter (iter=1000)
c* check the bisection condition
	if(f(a)*f(b).gt.0.0) then
           iflag = 0
           return
	end if
	i=0
	xl = a
	xr = b
	do while (abs(xr-xl).gt.eps)
           i=i+1
           x0=(xr+xl)/2
           if(f(xl)*f(x0).le.0.0) then
              xr = x0
              else
              xl=x0
           end if
           if(i.ge.iter) exit
	end do
	Root=x0
	iflag=1
	return
        end