function opt_dfunc(p)

   use nrtype
   use precision_def

   implicit none

    real(rp), dimension(:), intent(inout) :: p
    reAL(rp), DIMENSION(size(p)) :: opt_dfunc

    interface
      function opt_func(p)
      use nrtype
      use precision_def
      implicit none
      real(rp) opt_func
      real(rp), dimension(:), intent(in) :: p
      end function opt_func
    end interface

    real(rp), dimension(size(p)) ::  xx, g
    real(rp) fp, f0
    integer n
    integer :: way = 1, IGRAD = 18
    real(rp) :: scale = 1.e-5

    f0 = opt_func(p)

    xx = p

     way = -way

    do n=1, ubound(p, 1)

     p(n) = xx(n) + 0.5**igrad*way
     fp = opt_func(p)
     G(n) = way*(fp-f0)*2.**IGRAD
     p(n) = xx(n)

     end do
     
     if(maxval(G) > 1.e18 .or. minval(G) < -1.e18)G=G*1.e-10
     opt_dfunc = G /sqrt(dot_product(G,G)) * scale
   
end function