qe_erfc Function

  function qe_erfc(x)
    !---------------------------------------------------------------------
    !
    !! erfc(x) = 1-erf(x)  - See comments in erf
    !
    use w90_constants, only: dp
    implicit none
    real(kind=dp), intent(in) :: x
    real(kind=dp)            :: qe_erfc
    real(kind=dp) :: ax, x2, xm2, p2(8), q2(8), p3(5), q3(5), pim1
    data p2/3.004592610201616E2_dp, 4.519189537118719E2_dp, &
      3.393208167343437E2_dp, 1.529892850469404E2_dp, &
      4.316222722205674E1_dp, 7.211758250883094_dp, &
      5.641955174789740E-1_dp, -1.368648573827167E-7_dp/
    data q2/3.004592609569833E2_dp, 7.909509253278980E2_dp, &
      9.313540948506096E2_dp, 6.389802644656312E2_dp, &
      2.775854447439876E2_dp, 7.700015293522947E1_dp, &
      1.278272731962942E1_dp, 1.000000000000000_dp/
    data p3/-2.996107077035422E-3_dp, -4.947309106232507E-2_dp, &
      -2.269565935396869E-1_dp, -2.786613086096478E-1_dp, &
      -2.231924597341847E-2_dp/
    data q3/1.062092305284679E-2_dp, 1.913089261078298E-1_dp, &
      1.051675107067932_dp, 1.987332018171353_dp, &
      1.000000000000000_dp/

    data pim1/0.56418958354775629_dp/
    !        ( pim1= sqrt(1/pi) )
    ax = abs(x)
    if (ax > 26.0_dp) then
      !
      !  erfc(26.0)=10^(-296); erfc( 9.0)=10^(-37);
      !
      qe_erfc = 0.0_dp
    elseif (ax > 4.0_dp) then
      x2 = x**2
      xm2 = (1.0_dp/ax)**2
      qe_erfc = (1.0_dp/ax)*exp(-x2)*(pim1 + xm2*(p3(1) + xm2*(p3(2) + xm2*(p3(3) + xm2*(p3(4) + xm2*p3(5)))))/ &
                                      (q3(1) + xm2*(q3(2) + xm2*(q3(3) + xm2*(q3(4) + xm2*q3(5))))))
    elseif (ax > 0.47_dp) then
      x2 = x**2
      qe_erfc = exp(-x2)* &
                (p2(1) + ax*(p2(2) + ax*(p2(3) + ax*(p2(4) + ax*(p2(5) + ax*(p2(6) + ax*(p2(7) + ax*p2(8))))))))/ &
                (q2(1) + ax*(q2(2) + ax*(q2(3) + ax*(q2(4) + ax*(q2(5) + ax*(q2(6) + ax*(q2(7) + ax*q2(8))))))))
    else
      qe_erfc = 1.0_dp - qe_erf(ax)
    endif
    !
    ! erf(-x)=-erf(x)  =>  erfc(-x) = 2-erfc(x)
    !
    if (x < 0.0_dp) qe_erfc = 2.0_dp - qe_erfc
    !
    return
  end function qe_erfc