subroutine grabgr(r,gradup,graddn,i,imin,imax,done,gagd,gagup,
     j     gagdn)
c  calculates grad(abs(grad(density)))
c  input r            : radial coordinate (array)
c  input gradup       : grad(density), spin up (array)
c  input graddn       : grad(density), spin down (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 abs(grad(density)) is splined
c  output gagd        : grad(abs(grad(density)))
c  output gagup,gagdn : grad(abs(grad(density))), spin up and down
      implicit logical (a-z)
      logical done
      double precision r,gradup,graddn,gagd,gagup,gagdn,temp,gag2,deriv
      integer i,imin,imax,j
      common/ggagag/temp,gag2
      dimension r(561),gradup(561),graddn(561)
      dimension temp(2,561),gag2(2,561)
      if (.not.done) then
         do 10 j = imin,imax
            temp(1,j) = dabs(gradup(j))
            temp(2,j) = dabs(graddn(j))
 10      continue
         call spline(r,imin+1,imax,temp,gag2)
         done = .true.
      endif 
      gagup = deriv(r,temp,gag2,1,i,imin,imax)
      gagdn = deriv(r,temp,gag2,-1,i,imin,imax)
      gagd = gagup + gagdn
      return
      end