cerror Subroutine

   subroutine cerror(z,Cer)

      complex(wp),intent(in) :: z !! Complex argument
      complex(wp),intent(out) :: Cer !! `erf(z)`

      complex(wp) :: c0 , cl , cr , cs , z1
      integer :: k
      real(wp) :: a0

      a0 = abs(z)
      c0 = exp(-z*z)
      z1 = z
      if ( real(z,wp)<0.0_wp ) z1 = -z
      !
      !       Cutoff radius R = 4.36; determined by balancing rounding error
      !       and asymptotic expansion error, see below.
      !
      !       The resulting maximum global accuracy expected is around 1e-8
      !
      if ( a0<=4.36_wp ) then
         !
         !          Rounding error in the Taylor expansion is roughly
         !
         !          ~ R*R * EPSILON * R**(2 R**2) / (2 R**2 Gamma(R**2 + 1/2))
         !
         cs = z1
         cr = z1
         do k = 1 , 120
            cr = cr*z1*z1/(k+0.5_wp)
            cs = cs + cr
            if ( abs(cr/cs)<1.0e-15_wp ) exit
         enddo
         Cer = 2.0_wp*c0*cs/sqrtpi
      else
         cl = 1.0_wp/z1
         cr = cl
         !
         !          Asymptotic series; maximum K must be at most ~ R^2.
         !
         !          The maximum accuracy obtainable from this expansion is roughly
         !
         !          ~ Gamma(2R**2 + 2) / (
         !                   (2 R**2)**(R**2 + 1/2) Gamma(R**2 + 3/2) 2**(R**2 + 1/2))
         !
         do k = 1 , 20
            cr = -cr*(k-0.5_wp)/(z1*z1)
            cl = cl + cr
            if ( abs(cr/cl)<1.0e-15_wp ) exit
         enddo
         Cer = 1.0_wp - c0*cl/sqrtpi
      endif
      if ( real(z,wp)<0.0_wp ) Cer = -Cer

   end subroutine cerror