real function sqrt(x)
************************************************************************
*     (Square Root)
*     This function computes the square root of its argument via
*     a cubically-convergent iterative algorithm published by
*
*     Cleve Moler and Donald Morrison, "Replacing Square Roots by
*     Pythagorean Sums",  IBM J.  Research and  Development,  27,
*     577-581 (1983).
*
*     Augustin A. Dubrulle, "A Class of Numerical Methods for the
*     Computation of  Pythagorean  Sums",  IBM  J.  Research  and
*     Development, 27, 582-589 (1983).
*
*     The algorithm has been augmented by checks for negative, NaN,
*     and Infinity arguments.  These exceptional conditions generate
*     run-time NaN (negative or NaN x) or Infinity values that can be
*     trapped if necessary.
*
*     This algorithm is ABSYMAL for |x| << 1; for x = 0, after 10,000
*     iterations, it had converged only to p = 1.e-6.
*     (12-Feb-1994)
************************************************************************

      real ZERO
      parameter (ZERO = 0.0e+00)
      real ONE
      parameter (ONE = 1.0e+00)
      real TWO
      parameter (TWO = 2.0e+00)
      real FOUR
      parameter (FOUR = 4.0e+00)

      real p, r, s, x
      logical isNaN, isInf
      integer n

      if (x .lt. ZERO) then
           sqrt = ZERO
           sqrt = sqrt / sqrt
      else if (isNaN(x)) then
           sqrt = x + x
      else if (isInf(x)) then
           sqrt = x + x
      else

           n = 0

           p = ONE
           r = x - ONE
   10      if ((ONE + r) .ne. ONE) then
                s = r / (FOUR + r)
                p = p + TWO * s * p
                r = r * (s/(ONE + TWO * s))**2

                n = n + 1
                if (n .gt. 5) print *,'x =',x,' s = ',s,' p =',p,
     x               ' r =',r
                if (n .gt. 10) then
                     print *,
     x                   '---------- mmsqrt() not converging ----------'
                else
                     go to 10
                endif
           endif
           sqrt = p
      endif

      end