complex(kind(1.d0)) function zlnrel (z)

        implicit none
  c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
  c
  c clnrel(z) = clog(1+z) with relative error accuracy near z = 0.
  c let   rho = cabs(z)  and
  c       r**2 = cabs(1+z)**2 = (1+x)**2 + y**2 = 1 + 2*x + rho**2 .
  c now if rho is small we may evaluate clnrel(z) accurately by
  c       clog(1+z) = cmplx (alog(r), carg(1+z))
  c                 = cmplx (0.5*alog(r**2), carg(1+z))
  c                 = cmplx (0.5*alnrel(2*x+rho**2), carg(1+z))
  c

        complex(kind(1.d0)) z
        real(kind(1.d0)) dlnrel,d1mach,zarg,sqeps,x,rho

        external dlnrel, zarg, d1mach

        data sqeps /0.0d0/

  c
        if (sqeps.eq.0.) sqeps = sqrt (d1mach(4))
  c
        if ( abs(1.+z).lt.sqeps) call seteru (
       1  54hzlnrel  answer lt half precision because z too near -1, 54,
       2  1, 1)
  c
        rho = abs(z)
        if (rho.gt.0.375) zlnrel = log (1.0+z)
        if (rho.gt.0.375) return
  c
        x = real(z)

        zlnrel = cmplx (0.5*dlnrel(2.*x+rho**2), zarg(1.0+z), kind(1.d0))

  c
        return
        end