bracket Subroutine

subroutine bracket(ax,bx,cx,fa,fb,fc,func)
  real(8), intent(inout) :: ax,bx
  real(8), intent(out) :: cx,fa,fb,fc
  !...the first parameter is the default ratio by which successive intervals
  !   are magnified; the second is the maximum magnification allowed for a
  !   parabolic-fit step
  real(8), parameter :: gold=1.618034d0,glimit=100.d0,tiny=1.d-20
  real(8) :: fu,q,r,u,ulim
  interface
     function func(x)
       real(8) :: x
       real(8) :: func
     end function func
  end interface
  fa=func(ax)
  fb=func(bx)
  if (fb > fa) then
     call swap(ax,bx)
     call swap(fa,fb)
  end if
  cx=bx+GOLD*(bx-ax)
  fc=func(cx)
  do
     if (fb < fc) RETURN
     r=(bx-ax)*(fb-fc)
     q=(bx-cx)*(fb-fa)
     u=bx-((bx-cx)*q-(bx-ax)*r)/(2.0*sign(max(abs(q-r),TINY),q-r))
     ulim=bx+GLIMIT*(cx-bx)
     if ((bx-u)*(u-cx) > 0.d0) then
        fu=func(u)
        if (fu < fc) then
           ax=bx
           fa=fb
           bx=u
           fb=fu
           RETURN
        else if (fu > fb) then
           cx=u
           fc=fu
           RETURN
        end if
        u=cx+GOLD*(cx-bx)
        fu=func(u)
     else if ((cx-u)*(u-ulim) > 0.d0) then
        fu=func(u)
        if (fu < fc) then
           bx=cx
           cx=u
           u=cx+GOLD*(cx-bx)
           call shft(fb,fc,fu,func(u))
        end if
     else if ((u-ulim)*(ulim-cx) >= 0.d0) then
        u=ulim
        fu=func(u)
     else
        u=cx+GOLD*(cx-bx)
        fu=func(u)
     end if
     call shft(ax,bx,cx,u)
     call shft(fa,fb,fc,fu)
  end do
contains
  subroutine swap(a,b)
    real(8), intent(inout) :: a,b
    real(8) :: dum
    dum=a
    a=b
    b=dum
  end subroutine swap
  !-------------------
  subroutine shft(a,b,c,d)
    real(8), intent(out) :: a
    real(8), intent(inout) :: b,c
    real(8), intent(in) :: d
    a=b
    b=c
    c=d
  end subroutine shft
end subroutine bracket