!******************************************************************************
! Complementary error function
! From GNU scientific library 0.6, translated from c to fortran.
!******************************************************************************

      function erfc11(x) result(erfc)

      use precision

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      real(wp)              :: ax,t,num,den
      integer               :: i
      real(wp), parameter   :: p(0:5) = (/ &
                                      0.0000000297886562639399288862e+08_wp, &
                                      0.0000007409740605964741794425e+07_wp, &
                                      0.0000061602098531096305440906e+06_wp, &
                                      0.0005019049726784267463450058e+04_wp, &
                                      0.1275366644729965952479585264e+01_wp, &
                                      0.5641895835477550741253201704e+00_wp  /)
      real(wp), parameter   :: q(0:6) = (/ &
                                      0.0000000033690752069827527677e+09_wp, &
                                      0.0000009608965327192787870698e+07_wp, &
                                      0.0001708144074746600431571095e+05_wp, &
                                      0.0120489519278551290360340491e+03_wp, &
                                      0.9396034016235054150430579648e+01_wp, &
                                      0.2260528520767326969591866945e+01_wp, &
                                      0.1000000000000000000000000000e+01_wp  /)
      real(wp), parameter   :: c1(20) = (/ &
                                   +0.106073416421769980345174155056e+01_wp, &
                                   -0.042582445804381043569204735291e+01_wp, &
                                   +0.004955262679620434040357683080e+01_wp, &
                                   +0.000449293488768382749558001242e+01_wp, &
                                   -0.000129194104658496953494224761e+01_wp, &
                                   -0.000001836389292149396270416979e+01_wp, &
                                   +0.000002211114704099526291538556e+01_wp, &
                                   -0.523337485234257134673693179020e-06_wp, &
                                   -0.278184788833537885382530989578e-06_wp, &
                                   +0.141158092748813114560316684249e-07_wp, &
                                   +0.272571296330561699984539141865e-08_wp, &
                                   -0.206343904872070629406401492476e-09_wp, &
                                   -0.214273991996785367924201401812e-10_wp, &
                                   +0.222990255539358204580285098119e-11_wp, &
                                   +0.136250074650698280575807934155e-12_wp, &
                                   -0.195144010922293091898995913038e-13_wp, &
                                   -0.685627169231704599442806370690e-15_wp, &
                                   +0.144506492869699938239521607493e-15_wp, &
                                   +0.245935306460536488037576200030e-17_wp, &
                                   -0.929599561220523396007359328540e-18_wp  /)
      real(wp), parameter   :: c2(25) = (/ &
                                   +0.044045832024338111077637466616e+01_wp, &
                                   -0.143958836762168335790826895326e+00_wp, &
                                   +0.044786499817939267247056666937e+00_wp, &
                                   -0.013343124200271211203618353102e+00_wp, &
                                   +0.003824682739750469767692372556e+00_wp, &
                                   -0.001058699227195126547306482530e+00_wp, &
                                   +0.000283859419210073742736310108e+00_wp, &
                                   -0.000073906170662206760483959432e+00_wp, &
                                   +0.000018725312521489179015872934e+00_wp, &
                                   -0.462530981164919445131297264430e-05_wp, &
                                   +0.111558657244432857487884006422e-05_wp, &
                                   -0.263098662650834130067808832725e-06_wp, &
                                   +0.607462122724551777372119408710e-07_wp, &
                                   -0.137460865539865444777251011793e-07_wp, &
                                   +0.305157051905475145520096717210e-08_wp, &
                                   -0.665174789720310713757307724790e-09_wp, &
                                   +0.142483346273207784489792999706e-09_wp, &
                                   -0.300141127395323902092018744545e-10_wp, &
                                   +0.622171792645348091472914001250e-11_wp, &
                                   -0.126994639225668496876152836555e-11_wp, &
                                   +0.255385883033257575402681845385e-12_wp, &
                                   -0.506258237507038698392265499770e-13_wp, &
                                   +0.989705409478327321641264227110e-14_wp, &
                                   -0.190685978789192181051961024995e-14_wp, &
                                   +0.350826648032737849245113757340e-15_wp  /)
      real(wp), parameter :: c3(20) = (/ &
                                   +0.111684990123545698684297865808e+01_wp, &
                                   +0.003736240359381998520654927536e+00_wp, &
                                   -0.000916623948045470238763619870e+00_wp, &
                                   +0.000199094325044940833965078819e+00_wp, &
                                   -0.000040276384918650072591781859e+00_wp, &
                                   +0.776515264697061049477127605790e-05_wp, &
                                   -0.144464794206689070402099225301e-05_wp, &
                                   +0.261311930343463958393485241947e-06_wp, &
                                   -0.461833026634844152345304095560e-07_wp, &
                                   +0.800253111512943601598732144340e-08_wp, &
                                   -0.136291114862793031395712122089e-08_wp, &
                                   +0.228570483090160869607683087722e-09_wp, &
                                   -0.378022521563251805044056974560e-10_wp, &
                                   +0.617253683874528285729910462130e-11_wp, &
                                   -0.996019290955316888445830597430e-12_wp, &
                                   +0.158953143706980770269506726000e-12_wp, &
                                   -0.251045971047162509999527428316e-13_wp, &
                                   +0.392607828989125810013581287560e-14_wp, &
                                   -0.607970619384160374392535453420e-15_wp, &
                                   +0.912600607264794717315507477670e-16_wp  /)


      ax = abs(x)
      if (ax == 0.0_wp) then
          stop 'GSL_ERFC: zero argument'
      elseif (ax <= 1.0_wp) then
          t = 2.0_wp*ax - 1.0_wp
          erfc = evalcs(t,c1)
      elseif (ax <= 5.0_wp) then
          t = 0.5_wp*(ax-3.0_wp)
          erfc = evalcs(t,c2) * exp(-x**2)
      elseif (ax < 10.0_wp) then
          t = (2.0_wp*x - 15.0_wp)/5.0_wp
          erfc = evalcs(t,c3) * exp(-x**2)/x
      else
          num = p(5)
          do i=4,0,-1
              num = x*num + p(i)
          enddo
          den = q(6)
          do i=5,0,-1
              den = x*den + q(i)
          enddo
          erfc = num/den * exp(-x**2)
      endif

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

      contains

!*******************************************************************************

      function evalcs(x,cs)

      implicit none
      real(wp), intent(in)  :: x,cs(:)
      real(wp)              :: evalcs
      real(wp)              :: b0,b1,b2,twox
      integer               :: i


      b0 = 0.0_wp
      b1 = 0.0_wp
      twox = 2.0_wp*x
      do i=size(cs),1,-1
          b2 = b1
          b1 = b0
          b0 = twox*b1 - b2 + cs(i)
      enddo
      evalcs = (b0-b2)/2.0_wp

      end function evalcs

!*******************************************************************************

      end function erfc11
