!*******************************************************************************
! Complementary error function in double precision
! COPYRIGHT : M. Mori  AUG. 30 1990  V.1-2
!*******************************************************************************

      function erfc08(x) result(erfc)

      use precision

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      integer               :: i,na,nm
      real(wp)              :: v,xv,y
      real(wp), parameter   :: ph = 1.772453850905516_wp
      real(wp), parameter   :: cmc = 2.0_wp / ph
      real(wp), parameter   :: h = 0.5_wp
      real(wp), parameter   :: cxc = 2.0_wp * h / 3.141592653589793_wp
      real(wp), parameter   :: cxi = 4.0_wp / cxc
      real(wp), parameter   :: cm(0:5) = (/  0.1000000000000000e+01_wp,  &
                                            -0.3333333333333333e+00_wp,  &
                                             0.1000000000000000e+00_wp,  &
                                            -0.2380952380952381e-01_wp,  &
                                             0.4629629629629630e-02_wp,  &
                                            -0.7575757575757575e-03_wp  /)
      real(wp), parameter   :: cx(1:13) = (/ 0.7788007830714048e+00_wp,  &
                                             0.3678794411714423e+00_wp,  &
                                             0.1053992245618643e+00_wp,  &
                                             0.1831563888873418e-01_wp,  &
                                             0.1930454136227709e-02_wp,  &
                                             0.1234098040866796e-03_wp,  &
                                             0.4785117392129009e-05_wp,  &
                                             0.1125351747192591e-06_wp,  &
                                             0.1605228055185612e-08_wp,  &
                                             0.1388794386496402e-10_wp,  &
                                             0.7287724095819692e-13_wp,  &
                                             0.2319522830243569e-15_wp,  &
                                             0.4477732441718302e-18_wp  /)
       real(wp), parameter  :: cq(1:13) = (/  0.25_wp,  1.00_wp,  2.25_wp,  &
                                              4.00_wp,  6.25_wp,  9.00_wp,  &
                                             12.25_wp, 16.00_wp, 20.25_wp,  &
                                             25.00_wp, 30.25_wp, 36.00_wp,  &
                                             42.25_wp  /)
       real(wp), parameter  :: ca(0:5) = (/   1.0_wp,  -1.0_wp,    3.0_wp,  &
                                            -15.0_wp, 105.0_wp, -945.0_wp  /)



      if (x <= 0.0_wp) then
          xv = -x
      else
          xv = x
      endif

      if (xv <= 0.1_wp) then
          y = xv**2
          nm = ubound(cm,1)
          v = cm(nm)
          do i=nm-1,0,-1
              v = cm(i) + y*v
          enddo
          erfc = 1.0_wp - cmc*xv*v
      elseif (xv <= 100.0_wp) then
          y = xv**2
          v = 1.0_wp/(2.0_wp*y)
          do i=1,ubound(cx,1)
              v = v + cx(i) / (cq(i)+y)
          enddo
          v = cxc*xv*exp(-y)*v
          if (xv < 6.0_wp) v = v - 2.0_wp/(exp(cxi*xv) - 1.0_wp)
          erfc = v
      else
          y = 1.0_wp/(2.0_wp*xv*xv)
          na = ubound(ca,1)
          v = ca(na)
          do i=na-1,0,-1
             v = ca(i) + y*v
          enddo
          v = exp(-xv**2)/(ph*xv)*v
          erfc = v
      endif

      if (x < 0)  erfc = 2.0_wp - erfc

      end function erfc08