! Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).


      function erfc09(x) result(erfc)

      use precision

      implicit none
      real(wp), intent(in)     :: x
      real(wp)                 :: erfc,y
      real(wp), parameter      :: pv = 1.26974899965115684e+01_wp
      real(wp), parameter      :: ph = 6.10399733098688199e+00_wp
      real(wp), parameter      :: p0 = 2.96316885199227378e-01_wp
      real(wp), parameter      :: p1 = 1.81581125134637070e-01_wp
      real(wp), parameter      :: p2 = 6.81866451424939493e-02_wp
      real(wp), parameter      :: p3 = 1.56907543161966709e-02_wp
      real(wp), parameter      :: p4 = 2.21290116681517573e-03_wp
      real(wp), parameter      :: p5 = 1.91395813098742864e-04_wp
      real(wp), parameter      :: p6 = 9.71013284010551623e-06_wp
      real(wp), parameter      :: p7 = 1.66642447174307753e-07_wp
      real(wp), parameter      :: q0 = 6.12158644495538758e-02_wp
      real(wp), parameter      :: q1 = 5.50942780056002085e-01_wp
      real(wp), parameter      :: q2 = 1.53039662058770397e+00_wp
      real(wp), parameter      :: q3 = 2.99957952311300634e+00_wp
      real(wp), parameter      :: q4 = 4.95867777128246701e+00_wp
      real(wp), parameter      :: q5 = 7.41471251099335407e+00_wp
      real(wp), parameter      :: q6 = 1.04765104356545238e+01_wp
      real(wp), parameter      :: q7 = 1.48455557345597957e+01_wp


      y = x**2
      y = exp(-y)*x*(  p7/(y + q7)   &
                     + p6/(y + q6)   &
                     + p5/(y + q5)   &
                     + p4/(y + q4)   &
                     + p3/(y + q3)   &
                     + p2/(y + q2)   &
                     + p1/(y + q1)   &
                     + p0/(y + q0) )
      if (x < ph) y = y + 2.0_wp/(exp(pv*x) + 1.0_wp)
      erfc = y

      end function erfc09