function erf (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension erfcs(13)
      external alog, csevl, erfc, inits, r1mach, sqrt
c
c series for erf        on the interval  0.          to  1.00000d+00
c                                        with weighted error   7.10e-18
c                                         log weighted error  17.15
c                               significant figures required  16.31
c                                    decimal places required  17.71
c
      data erf cs( 1) /   -.0490461212 34691808e0 /
      data erf cs( 2) /   -.1422612051 0371364e0 /
      data erf cs( 3) /    .0100355821 87599796e0 /
      data erf cs( 4) /   -.0005768764 69976748e0 /
      data erf cs( 5) /    .0000274199 31252196e0 /
      data erf cs( 6) /   -.0000011043 17550734e0 /
      data erf cs( 7) /    .0000000384 88755420e0 /
      data erf cs( 8) /   -.0000000011 80858253e0 /
      data erf cs( 9) /    .0000000000 32334215e0 /
      data erf cs(10) /   -.0000000000 00799101e0 /
      data erf cs(11) /    .0000000000 00017990e0 /
      data erf cs(12) /   -.0000000000 00000371e0 /
      data erf cs(13) /    .0000000000 00000007e0 /
c
      data sqrtpi /1.772453850 9055160e0/
      data nterf, xbig, sqeps / 0, 0., 0./
c
      if (nterf.ne.0) go to 10
      nterf = inits (erfcs, 13, 0.1*r1mach(3))
      xbig = sqrt(-alog(sqrtpi*r1mach(3)))
      sqeps = sqrt(2.0*r1mach(3))
c
 10   y = abs(x)
      if (y.gt.1.) go to 20
c
c erf(x) = 1. - erfc(x) for -1. .le. x .le. 1.
c
      if (y.le.sqeps) erf = 2.0*x/sqrtpi
      if (y.gt.sqeps) erf = x*(1.0 + csevl(2.*x**2-1., erfcs, nterf))
      return
c
c erf(x) = 1. - erfc(x) for  abs(x) .gt. 1.
c
 20   if (y.le.xbig) erf = sign (1.0-erfc(y), x)
      if (y.gt.xbig) erf = sign (1.0, x)
c
      return
      end