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

      end