function sqrt (x)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension sqrt2(3)
      external alog, r1mach, r9pak
      data sqrt2(1) / 0.7071067811 8654752e0 /
      data sqrt2(2) / 1.0 e0 /
      data sqrt2(3) / 1.4142135623 7309505 e0 /
c
      data niter / 0 /
c
      if (niter.eq.0) niter = 1.443*alog(-0.104*alog(0.1*r1mach(3)))+ 1.
c
      if (x.le.0.) go to 20
c
      call r9upak (x, y, n)
      ixpnt = n/2
      irem = n - 2*ixpnt + 2
c
c the approximation below has accuracy of 4.16 digits.
      sqrt = .261599e0 + y*(1.114292e0 + y*(-.516888e0 + y*.141067e0))
c
      do 10 iter=1,niter
        sqrt = sqrt + 0.5*(y - sqrt**2) / sqrt
 10   continue
c
      sqrt = r9pak (sqrt2(irem)*sqrt, ixpnt)
      return
c
 20   if (x.lt.0.) call seteru (21hsqrt    x is negative, 21, 1, 1)
      sqrt = 0.0
      return
c
      end