cerf Subroutine

   subroutine cerf(z,Cer,Cder)

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

      real(wp) :: c0 , cs , ei1 , ei2 , er , er0 , er1 , &
                  er2 , eri , err , r , ss , w , w1 , w2 , x ,&
                  x2 , y
      integer :: k , n

      real(wp),parameter :: eps = 1.0e-12_wp

      x = real(z,wp)
      y = aimag(z)
      x2 = x*x
      if ( x<=3.5_wp ) then
         er = 1.0_wp
         r = 1.0_wp
         w = 0.0_wp
         do k = 1 , 100
            r = r*x2/(k+0.5_wp)
            er = er + r
            if ( abs(er-w)<=eps*abs(er) ) exit
            w = er
         enddo
         c0 = 2.0_wp/sqrtpi*x*exp(-x2)
         er0 = c0*er
      else
         er = 1.0_wp
         r = 1.0_wp
         do k = 1 , 12
            r = -r*(k-0.5_wp)/x2
            er = er + r
         enddo
         c0 = exp(-x2)/(x*sqrtpi)
         er0 = 1.0_wp - c0*er
      endif
      if ( y==0.0_wp ) then
         err = er0
         eri = 0.0_wp
      else
         cs = cos(2.0_wp*x*y)
         ss = sin(2.0_wp*x*y)
         er1 = exp(-x2)*(1.0_wp-cs)/(twopi*x)
         ei1 = exp(-x2)*ss/(twopi*x)
         er2 = 0.0_wp
         w1 = 0.0_wp
         do n = 1 , 100
            er2 = er2 + exp(-0.25_wp*n*n)/(n*n+4.0_wp*x2) &
                  *(2.0_wp*x-2.0_wp*x*cosh(n*y)*cs+n*sinh(n*y)*ss)
            if ( abs((er2-w1)/er2)<eps ) exit
            w1 = er2
         enddo
         c0 = 2.0_wp*exp(-x2)/pi
         err = er0 + er1 + c0*er2
         ei2 = 0.0_wp
         w2 = 0.0_wp
         do n = 1 , 100
            ei2 = ei2 + exp(-0.25_wp*n*n)/(n*n+4.0_wp*x2) &
                  *(2.0_wp*x*cosh(n*y)*ss+n*sinh(n*y)*cs)
            if ( abs((ei2-w2)/ei2)<eps ) exit
            w2 = ei2
         enddo
         eri = ei1 + c0*ei2
      endif
      Cer = cmplx(err,eri,kind=wp)
      Cder = 2.0_wp/sqrtpi*exp(-z*z)

   end subroutine cerf