!*******************************************************************************
! Complementary Error Function
!
!   erfc(z)=f(z)+(2*h/pi)*exp(-z^2)*z
!           *(exp(-(1*h/2)^2)/(z^2+(1*h/2)^2)
!           +exp(-(3*h/2)^2)/(z^2+(3*h/2)^2))
!           +exp(-(5*h/2)^2)/(z^2+(5*h/2)^2)
!           +...)
!       f(z)=2/(1+exp(2*pi*z/h)) , Re(z)+abs(Im(z))<pi/h
!       f(z)=0                   , Re(z)+abs(Im(z))>=pi/h
!
!   erfc(z)=g(z)+(2*h/pi)*exp(-z^2)*z
!           *(1/(2*z^2)
!           +exp(-(1*h)^2)/(z^2+(1*h)^2)
!           +exp(-(2*h)^2)/(z^2+(2*h)^2)
!           +exp(-(3*h)^2)/(z^2+(3*h)^2)
!           +...)
!       g(z)=2/(1-exp(2*pi*z/h)) , Re(z)+abs(Im(z))<pi/h
!       g(z)=0                   , Re(z)+abs(Im(z))>=pi/h
!
!   M. Mori, A Method for Evaluation of the Error Function
!   of Real and Complex Variable with High Relative Accuracy,
!   Publ. RIMS, Kyoto Univ. Vol. 19, 1983.
!*******************************************************************************

      function erfc06(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)     :: x
      real(wp)                 :: erfc
      real(wp)                 :: t,u
      real(wp), parameter      :: pa  =  3.97886080735226000e+00_wp
      real(wp), parameter      :: p00 =  2.75374741597376782e-01_wp
      real(wp), parameter      :: p01 =  4.90165080585318424e-01_wp
      real(wp), parameter      :: p02 =  7.74368199119538609e-01_wp
      real(wp), parameter      :: p03 =  1.07925515155856677e+00_wp
      real(wp), parameter      :: p04 =  1.31314653831023098e+00_wp
      real(wp), parameter      :: p05 =  1.37040217682338167e+00_wp
      real(wp), parameter      :: p06 =  1.18902982909273333e+00_wp
      real(wp), parameter      :: p07 =  8.05276408752910567e-01_wp
      real(wp), parameter      :: p08 =  3.57524274449531043e-01_wp
      real(wp), parameter      :: p09 =  1.66207924969367356e-02_wp
      real(wp), parameter      :: p10 = -1.19463959964325415e-01_wp
      real(wp), parameter      :: p11 = -8.38864557023001992e-02_wp
      real(wp), parameter      :: p12 =  2.49367200053503304e-03_wp
      real(wp), parameter      :: p13 =  3.90976845588484035e-02_wp
      real(wp), parameter      :: p14 =  1.61315329733252248e-02_wp
      real(wp), parameter      :: p15 = -1.33823644533460069e-02_wp
      real(wp), parameter      :: p16 = -1.27223813782122755e-02_wp
      real(wp), parameter      :: p17 =  3.83335126264887303e-03_wp
      real(wp), parameter      :: p18 =  7.73672528313526668e-03_wp
      real(wp), parameter      :: p19 = -8.70779635317295828e-04_wp
      real(wp), parameter      :: p20 = -3.96385097360513500e-03_wp
      real(wp), parameter      :: p21 =  1.19314022838340944e-04_wp
      real(wp), parameter      :: p22 =  1.27109764952614092e-03_wp


      t = pa/(pa + abs(x))
      u = t - 0.5_wp
      erfc = ((((((((((((((((((((((p22 * u + p21) * u + p20)   &
                                       * u + p19) * u + p18)   &
                                       * u + p17) * u + p16)   &
                                       * u + p15) * u + p14)   &
                                       * u + p13) * u + p12)   &
                                       * u + p11) * u + p10)   &
                                       * u + p09) * u + p08)   &
                                       * u + p07) * u + p06)   &
                                       * u + p05) * u + p04)   &
                                       * u + p03) * u + p02)   &
                                       * u + p01) * u + p00) * t * exp(-x**2)
      if (x < 0.0_wp) erfc = 2.0_wp - erfc

      end function erfc06