subroutine grad(r,db,i,imin,imax,done,gradd,gradup,graddn)
c  grad calculates the derivative of db with respect to r, at
c  position r(i).
c  input r              : radial coordinate (array)
c  input db             : spin-up and down densities (array)
c  input i              : counter, i.e. current position = r(i)
c  input imin           : min value of i
c  input imax           : max value of i
c  in/output done       : true if density is splined
c  output gradd         : grad(density)
c  output gradup,graddn : grad(density), spin up and down
      implicit logical (a-z)
      logical done
      double precision r,db,gradd,gradup,graddn,db2,deriv
      integer i,imin,imax
      dimension r(561),db(2,561),db2(2,561)
      common/ggagra/db2
      if (.not.done) then
         call spline(r,imin+1,imax,db,db2)
         done = .true.
      endif 
      gradup = deriv(r,db,db2,1,i,imin,imax)
      graddn = deriv(r,db,db2,-1,i,imin,imax)
      gradd = gradup + graddn
      return
      end