double precision function dnxtaft(x,y)
************************************************************************
*     Return the next floating-point number after x in the
*     direction of y.  If y = x, then dnxtaft(x,x) returns
*     x + epsilon.
*     (13-Nov-1990)
************************************************************************
      double precision x,y,temp,eps
      if (x .gt. y) then
         eps = x
 10      temp = x - eps/2.0d+00
         if (temp .eq. x) goto 20
         eps = eps/2.0d+00
         goto 10
 20      dnxtaft = x - eps
      else
         eps = x
 30      temp = x + eps/2.0d+00
         if (temp .eq. x) goto 40
         eps = eps/2.0d+00
         goto 30
 40      dnxtaft = x + eps
      endif

      end