!*******************************************************************************
! 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 erfc07(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)     :: x
      real(wp)                 :: erfc
      real(wp)                 :: t,u,v
      integer                  :: i
      real(wp), parameter      :: pa  =  3.97886080735226000e+00_wp
      real(wp), parameter      :: p(0:22) = (/ 2.75374741597376782e-01_wp,  &
                                               4.90165080585318424e-01_wp,  &
                                               7.74368199119538609e-01_wp,  &
                                               1.07925515155856677e+00_wp,  &
                                               1.31314653831023098e+00_wp,  &
                                               1.37040217682338167e+00_wp,  &
                                               1.18902982909273333e+00_wp,  &
                                               8.05276408752910567e-01_wp,  &
                                               3.57524274449531043e-01_wp,  &
                                               1.66207924969367356e-02_wp,  &
                                              -1.19463959964325415e-01_wp,  &
                                              -8.38864557023001992e-02_wp,  &
                                               2.49367200053503304e-03_wp,  &
                                               3.90976845588484035e-02_wp,  &
                                               1.61315329733252248e-02_wp,  &
                                              -1.33823644533460069e-02_wp,  &
                                              -1.27223813782122755e-02_wp,  &
                                               3.83335126264887303e-03_wp,  &
                                               7.73672528313526668e-03_wp,  &
                                              -8.70779635317295828e-04_wp,  &
                                              -3.96385097360513500e-03_wp,  &
                                               1.19314022838340944e-04_wp,  &
                                               1.27109764952614092e-03_wp   /)

      t = pa/(pa + abs(x))
      u = t - 0.5_wp
      v = p(22)
      do i=21,0,-1
          v = u*v + p(i)
      enddo
      erfc = v*t*exp(-x**2)
      if (x < 0.0_wp) erfc = 2.0_wp - erfc

      end function erfc07