dbrent_nograd Subroutine

subroutine dbrent_nograd(func,xmin,brack,tol,niter)
  interface
     function func(x)
       real(8) :: x
       real(8) :: func
     end function func
  end interface
  real(8),intent(inout)         :: xmin
  real(8),dimension(:),optional :: brack
  real(8),optional              :: tol
  integer,optional              :: niter
  real(8)                       :: tol_
  integer                       :: niter_
  integer                       :: iter
  real(8)                       :: ax,xx,bx,fa,fx,fb,fret
  !
  tol_=1d-9;if(present(tol))tol_=tol
  Niter_=200;if(present(Niter))Niter_=Niter
  !
  if(present(brack))then
     select case(size(brack))
     case(1)
        stop "Brent error: calling brent with size(brack)==1. None or two points are necessary."
     case(2)
        ax = brack(1)
        xx = brack(2)
        call bracket(ax,xx,bx,fa,fx,fb,func)
     case (3)
        ax = brack(1)
        xx = brack(2)
        bx = brack(3)
     end select
  else
     ax=0d0
     xx=1d0
     call bracket(ax,xx,bx,fa,fx,fb,func)
  endif
  fret=dbrent_optimize(ax,xx,bx,func,dfunc,tol_,niter_,xmin)
contains
  function dfunc(x)
    real(8) :: x
    real(8) :: dfunc
    call fgradient_func(func,x,dfunc)
  end function dfunc
  !
  subroutine fgradient_func(funcv,x,fjac,epsfcn)
    implicit none
    interface 
       function funcv(x)
         real(8) :: x
         real(8) :: funcv
       end function funcv
    end interface
    integer          ::  n
    real(8),intent(in) ::  x
    real(8)            ::  x_
    real(8)          ::  fvec
    real(8)          ::  fjac
    real(8),optional ::  epsfcn
    real(8)          ::  eps,eps_
    real(8)          ::  epsmch
    real(8)          ::  h,temp
    real(8)          ::  wa1
    real(8)          ::  wa2
    integer          :: i,j,k
    x_ = x
    eps_= 0.d0; if(present(epsfcn))eps_=epsfcn
    epsmch = epsilon(epsmch)
    eps  = sqrt(max(eps_,epsmch))
    !  Evaluate the function
    fvec = funcv(x_)
    temp = x_
    h    = eps*abs(temp)
    if(h==0d0) h = eps
    x_   = temp + h
    wa1  = funcv(x_)
    x_   = temp
    fjac = (wa1 - fvec)/h
  end subroutine fgradient_func
end subroutine dbrent_nograd