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

      function erfc02(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      integer               :: j
      real(wp)              :: a,b,y,xp(0:7)
      real(wp), parameter   :: pi = 3.1415926535897932384626433832795028_wp
      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  /)


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

      end function erfc02