!*******************************************************************************
! Complementary error-function  erfc(x) = 1 - erf(x)
! Reference: W.J. Cody, Mathematics of Computation 22 (1969), 631-637
! Taken from CERN library.
!*******************************************************************************

      function erfc01(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      integer               :: j
      real(wp)              :: a,b,v,y
      real(wp), parameter   :: isqrtpi = 0.56418958354775629_wp   ! 1/sqrt(pi)
      real(wp), parameter   :: p1(0:3) = (/  2.426679552305318e+2_wp, &
                                             2.197926161829415e+1_wp, &
                                             6.996383488619136e+0_wp, &
                                            -3.560984370181538e-2_wp  /)
      real(wp), parameter   :: q1(0:3) = (/  2.150588758698612e+2_wp, &
                                             9.116490540451490e+1_wp, &
                                             1.508279763040779e+1_wp, &
                                             1.000000000000000e+0_wp  /)
      real(wp), parameter   :: p2(0:7) = (/  3.004592610201616e+2_wp, &
                                             4.519189537118729e+2_wp, &
                                             3.393208167343437e+2_wp, &
                                             1.529892850469404e+2_wp, &
                                             4.316222722205674e+1_wp, &
                                             7.211758250883094e+0_wp, &
                                             5.641955174789740e-1_wp, &
                                            -1.368648573827167e-7_wp  /)
      real(wp), parameter   :: q2(0:7) = (/  3.004592609569833e+2_wp, &
                                             7.909509253278980e+2_wp, &
                                             9.313540948506096e+2_wp, &
                                             6.389802644656312e+2_wp, &
                                             2.775854447439876e+2_wp, &
                                             7.700015293522947e+1_wp, &
                                             1.278272731962942e+1_wp, &
                                             1.000000000000000e+0_wp  /)
      real(wp), parameter   :: p3(0:4) = (/ -2.996107077035422e-3_wp, &
                                            -4.947309106232507e-2_wp, &
                                            -2.269565935396869e-1_wp, &
                                            -2.786613086096478e-1_wp, &
                                            -2.231924597341847e-2_wp  /)
      real(wp), parameter   :: q3(0:4) = (/  1.062092305284679e-2_wp, &
                                             1.913089261078298e-1_wp, &
                                             1.051675107067932e+0_wp, &
                                             1.987332018171353e+0_wp, &
                                             1.000000000000000e+0_wp  /)


      v = abs(x)
      if (v <= 0.46875_wp) then
          y = v**2
          a = p1(3)
          b = q1(3)
          do j=2,0,-1
              a = a*y + p1(j)
              b = b*y + q1(j)
          enddo
          erfc = 1.0_wp - x*a/b
          return
      elseif (v <= 4.0_wp) then
          a = p2(7)
          b = q2(7)
          do j=6,0,-1
              a = a*v + p2(j)
              b = b*v + q2(j)
          enddo
          erfc = exp(-v**2)*a/b
      elseif (v <= 10.0_wp) then
          y = 1.0_wp/(v**2)
          a = p3(4)
          b = q3(4)
          do j=3,0,-1
              a = a*y + p3(j)
              b = b*y + q3(j)
          enddo
          erfc = exp(-v**2)*(isqrtpi+y*a/b)/v
      else
          erfc = 0.0_wp
      endif
      if (x <= 0.0_wp) erfc = 2.0_wp - erfc

      end function erfc01