c===========================================================
c     tbisect: Illustrates root finding using bisection 
c     routine 'bisect'.
c
c     Initial bracketing interval must be specified via the 
c     command-line, along with optional convergence criteria 
c     and output option.
c
c     This program also illustrates the general Fortran
c     techniques (briefly discussed previously) for: 
c
c     (1) Writing and using routines which take other routines
c     as arguments.
c     (2) Using a COMMON block to communicate information to
c     a routine in cases where the information cannot be 
c     passed via the argument list.
c     (3) Using an "INCLUDE" file (in this case 'comf.inc')
c     to ensure that the same common block structure is defined 
c     in all program units.
c
c     Currently set up for computing square roots i.e. 
c     solves 
c
c         f(x; a) = x**2 - a = 0
c
c     for 'a' specifed on command-line
c
c     Outputs a, approximate root (x*) and f(x*; a) on stdout.
c===========================================================
      program         tbisect
      implicit        none
c-----------------------------------------------------------
c     Declaration of the bisection routine.
c-----------------------------------------------------------
      real*8          bisect
c-----------------------------------------------------------
c     Name of the specific function whose root we seek.
c     Note use of 'external' to let compiler know 'fsqr'
c     is the name of a function, not a variable.
c-----------------------------------------------------------
      real*8          fsqr
      external        fsqr

      integer         i4arg,         iargc
      real*8          r8arg
c-----------------------------------------------------------
c     For use in detecting bad real*8 command-line value.
c-----------------------------------------------------------
      real*8          r8_never
      parameter     ( r8_never = -1.0d-60 )
c-----------------------------------------------------------
c     Use a common block to pass number whose square root 
c     is sought to external function 'fsqr'.
c-----------------------------------------------------------
      include        'comf.inc'
c-----------------------------------------------------------
c     Initial bracket, convergence tolerance and output
c     option from command-line; default value for conv.
c     tolerance.
c-----------------------------------------------------------
      real*8          xmin,           xmax,         xtol
      logical         trace

      real*8          default_xtol
      parameter     ( default_xtol = 1.0d-8 )
c-----------------------------------------------------------
c     Root and return code from bisection routine.
c-----------------------------------------------------------
      real*8          root
      integer         rc
c-----------------------------------------------------------
c     Argument parsing.
c-----------------------------------------------------------
      if( iargc() .lt. 3 ) go to 900
      a     = r8arg(1,r8_never)
      xmin  = r8arg(2,r8_never)
      xmax  = r8arg(3,r8_never)
      if( a .eq. r8_never  .or.  xmin .eq. r8_never  .or.  
     &    xmax .eq. r8_never ) go to 900

      xtol  = r8arg(4,default_xtol)
      trace = iargc() .gt. 4
c-----------------------------------------------------------
c     Invoke root finder then write a, sqrt(a), and residual
c     to standard output.
c-----------------------------------------------------------
      root = bisect(fsqr,xmin,xmax,xtol,trace,rc)
      if( rc .eq. 0 ) then
         write(*,*) a, root, fsqr(root)
      else
         write(0,*) 'tbisect: Bisection failed.'
      end if
c-----------------------------------------------------------
c     Normal exit.
c-----------------------------------------------------------
      stop

c-----------------------------------------------------------
c     Usage exit.
c-----------------------------------------------------------
 900  continue
         write(0,*) 'usage: tbisect <a> <xmin> <xmax> '//
     &              '[<xtol> <trace>]'
      stop

      end

c===========================================================
c     Function whose root is sought.  Again, note use of 
c     COMMON block to pass additional information (in this 
c     case 'a') to the routine.
c===========================================================
      real*8 function fsqr(x)
         implicit      none

         real*8        x

         include      'comf.inc'

         fsqr = x**2 - a

         return
      end