!*******************************************************************************
! Complementary Error Function
! Naval Surface Warfare Center, Dahlgren, Virginia, April 1992
!*******************************************************************************

      function erfc10(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      real(wp)              :: ax,t,w
      integer               :: i
      real(wp), parameter   :: a(21) = (/                                    &
                                  0.1283791670955125738961589031215e+00_wp,  &
                                 -0.3761263890318375246320529677070e+00_wp,  &
                                  0.1128379167095512573896158902931e+00_wp,  &
                                 -0.2686617064513125175943235372542e-01_wp,  &
                                  0.5223977625442187842111812447877e-02_wp,  &
                                 -0.8548327023450852832540164081187e-03_wp,  &
                                  0.1205533298178966425020717182498e-03_wp,  &
                                 -0.1492565035840625090430728526820e-04_wp,  &
                                  0.1646211436588924261080723578109e-05_wp,  &
                                 -0.1636584469123468757408968429674e-06_wp,  &
                                  0.1480719281587021715400818627811e-07_wp,  &
                                 -0.1229055530145120140800510155331e-08_wp,  &
                                  0.9422759058437197017313055084212e-10_wp,  &
                                 -0.6711366740969385085896257227159e-11_wp,  &
                                  0.4463222608295664017461758843550e-12_wp,  &
                                 -0.2783497395542995487275065856998e-13_wp,  &
                                  0.1634095572365337143933023780777e-14_wp,  &
                                 -0.9052845786901123985710019387938e-16_wp,  &
                                  0.4708274559689744439341671426731e-17_wp,  &
                                 -0.2187159356685015949749948252160e-18_wp,  &
                                  0.7043407712019701609635599701333e-20_wp   /)


      ax = abs(x)

!     abs(x) <= 1
      if (ax <= 1.0_wp) then
          t = x**2
          w = a(21)
          do i=20,1,-1
              w = t*w + a(i)
          enddo
          erfc = 0.5_wp + (0.5_wp - x*(1.0_wp + w))
          return
      endif

!     x < -1
      if (x <= 0.0_wp) then
          erfc = 2.0_wp
          if (x < -8.3_wp) return
          t = x**2
          erfc = 2.0_wp - exp(-t)*erfc0(ax)
          return
      endif

!     x > 1
      erfc = 0.0_wp
      if (x > 100.0_wp) return
      t = x**2
      if (t > -exparg(1)) return
      erfc = exp(-t) * erfc0(x)

      contains

!*******************************************************************************
! Evaluation of exp(x**2)*erfc(x) for x >= 1
! Written by Alfred H. Morris, Jr.
! Naval Surface Warfare Center, Dahlgren, Virginia, April 1992
!*******************************************************************************

      function erfc0(x)

      use precision, only: wp

      implicit none
      real(wp), intent(in) :: x
      real(wp)             :: erfc0
      real(wp)             :: t,u,v,z
!                             rpinv = 1/sqrt(pi)
      real(wp), parameter  :: rpinv = 0.56418958354775628694807945156077259_wp
      real(wp), parameter  :: p0 =  0.16506148041280876191828601e-03_wp,  &
                              p1 =  0.15471455377139313353998665e-03_wp,  &
                              p2 =  0.44852548090298868465196794e-04_wp,  &
                              p3 = -0.49177280017226285450486205e-05_wp,  &
                              p4 = -0.69353602078656412367801676e-05_wp,  &
                              p5 = -0.20508667787746282746857743e-05_wp,  &
                              p6 = -0.28982842617824971177267380e-06_wp,  &
                              p7 = -0.17272433544836633301127174e-07_wp
      real(wp), parameter  :: q0 =  1.00000000000000000000000000e+00_wp,  &
                              q1 =  0.16272656776533322859856317e+01_wp,  &
                              q2 =  0.12040996037066026106794322e+01_wp,  &
                              q3 =  0.52400246352158386907601472e+00_wp,  &
                              q4 =  0.14497345252798672362384241e+00_wp,  &
                              q5 =  0.25592517111042546492590736e-01_wp,  &
                              q6 =  0.26869088293991371028123158e-02_wp,  &
                              q7 =  0.13133767840925681614496481e-03_wp
      real(wp), parameter  :: r0 =  0.145589721275038539045668824025e+00_wp,  &
                              r1 = -0.273421931495426482902320421863e+00_wp,  &
                              r2 =  0.226008066916621506788789064272e+00_wp,  &
                              r3 = -0.163571895523923805648814425592e+00_wp,  &
                              r4 =  0.102604312032193978662297299832e+00_wp,  &
                              r5 = -0.548023266949835519254211506880e-01_wp,  &
                              r6 =  0.241432239725390106956523668160e-01_wp,  &
                              r7 = -0.822062115403915116036874169600e-02_wp,  &
                              r8 =  0.180296241564687154310619200000e-02_wp
      real(wp), parameter  :: a0 = -0.45894433406309678202825375e-03_wp,  &
                              a1 = -0.12281298722544724287816236e-01_wp,  &
                              a2 = -0.91144359512342900801764781e-01_wp,  &
                              a3 = -0.28412489223839285652511367e-01_wp,  &
                              a4 =  0.14083827189977123530129812e+01_wp,  &
                              a5 =  0.11532175281537044570477189e+01_wp,  &
                              a6 = -0.72170903389442152112483632e+01_wp,  &
                              a7 = -0.19685597805218214001309225e+01_wp,  &
                              a8 =  0.93846891504541841150916038e+01_wp
      real(wp), parameter  :: b0 =  1.00000000000000000000000000e+00_wp,  &
                              b1 =  0.25136329960926527692263725e+02_wp,  &
                              b2 =  0.15349442087145759184067981e+03_wp,  &
                              b3 = -0.29971215958498680905476402e+03_wp,  &
                              b4 = -0.33876477506888115226730368e+04_wp,  &
                              b5 =  0.28301829314924804988873701e+04_wp,  &
                              b6 =  0.22979620942196507068034887e+05_wp,  &
                              b7 = -0.24280681522998071562462041e+05_wp,  &
                              b8 = -0.36680620673264731899504580e+05_wp,  &
                              b9 =  0.42278731622295627627042436e+05_wp,  &
                              b10=  0.28834257644413614344549790e+03_wp,  &
                              b11=  0.70226293775648358646587341e+03_wp
      real(wp), parameter  :: c0 = -0.7040906288250128001000086e-04_wp,  &
                              c1 = -0.3858822461760510359506941e-02_wp,  &
                              c2 = -0.7708202127512212359395078e-01_wp,  &
                              c3 = -0.6713655014557429480440263e+00_wp,  &
                              c4 = -0.2081992124162995545731882e+01_wp,  &
                              c5 =  0.2898831421475282558867888e+01_wp,  &
                              c6 =  0.2199509380600429331650192e+02_wp,  &
                              c7 =  0.2907064664404115316722996e+01_wp,  &
                              c8 = -0.4766208741588182425380950e+02_wp
      real(wp), parameter  :: d0 =  1.0000000000000000000000000e+00_wp,  &
                              d1 =  0.5238852785508439144747174e+02_wp,  &
                              d2 =  0.9646843357714742409535148e+03_wp,  &
                              d3 =  0.7007152775135939601804416e+04_wp,  &
                              d4 =  0.8515386792259821780601162e+04_wp,  &
                              d5 = -0.1002360095177164564992134e+06_wp,  &
                              d6 = -0.2065250031331232815791912e+06_wp,  &
                              d7 =  0.5695324805290370358175984e+06_wp,  &
                              d8 =  0.6589752493461331195697873e+06_wp,  &
                              d9 = -0.1192930193156561957631462e+07_wp
      real(wp), parameter  :: e0 =  0.540464821348814822409610122136e+00_wp, &
                              e1 = -0.261515522487415653487049835220e-01_wp, &
                              e2 = -0.288573438386338758794591212600e-02_wp, &
                              e3 = -0.529353396945788057720258856000e-03_wp
      real(wp), parameter  :: s0  =  1.00000000000000000000e+00_wp,  &
                              s1  = -0.50000000000000000000e+00_wp,  &
                              s2  =  0.75000000000000000000e+00_wp,  &
                              s3  = -0.18750000000000000000e+01_wp,  &
                              s4  =  0.65625000000000000000e+01_wp,  &
                              s5  = -0.29531250000000000000e+02_wp,  &
                              s6  =  0.16242187500000000000e+03_wp,  &
                              s7  = -0.10557421875000000000e+04_wp,  &
                              s8  =  0.79180664062500000000e+04_wp,  &
                              s9  = -0.67303564453125000000e+05_wp,  &
                              s10 =  0.63938386230468750000e+06_wp,  &
                              s11 = -0.67135305541992187500e+07_wp,  &
                              s12 =  0.77205601373291015625e+08_wp


      if (x < 1.0_wp) then
          stop 'ERFC0: argument x < 1'
      elseif (x <= 2.0_wp) then
          u = ((((((p7*x + p6)*x + p5)*x + p4)*x + p3)*x + p2)*x + p1)*x + p0
          v = ((((((q7*x + q6)*x + q5)*x + q4)*x + q3)*x + q2)*x + q1)*x + q0
          t = (x-3.75_wp) / (x+3.75_wp)
          erfc0 = (((((((((u/v)*t + r8)*t + r7)*t + r6)*t + r5)  &
                               *t + r4)*t + r3)*t + r2)*t + r1)*t + r0
      elseif (x <= 4.0_wp) then
          z = 1.0_wp/(2.5_wp + x*x)
          u = (((((((a8*z + a7)*z + a6)*z + a5)*z + a4)  &
                       *z + a3)*z + a2)*z + a1)*z + a0
          v = ((((((((((b11*z + b10)*z + b9)*z + b8)*z + b7)*z + b6)  &
                           *z +  b5)*z + b4)*z + b3)*z + b2)*z + b1)*z + b0
          t = 13.0_wp*z - 1.0_wp
          erfc0 = ((((u/v)*t + e2)*t + e1)*t + e0) / x
      elseif (x < 50.0_wp) then
          z = 1.0_wp/(2.5_wp + x*x)
          u = (((((((c8*z + c7)*z + c6)*z + c5)*z + c4)  &
                       *z + c3)*z + c2)*z + c1)*z + c0
          v = ((((((((d9*z + d8)*z + d7)*z + d6)*z + d5)  &
                        *z + d4)*z + d3)*z + d2)*z + d1)*z + d0
          t = 13.0_wp*z - 1.0_wp
          erfc0 = (((((u/v)*t + e3)*t + e2)*t + e1)*t + e0) / x
      else
          t = 1.0_wp/(x*x)
          z = (((((((((((s12*t + s11)*t + s10)*t + s9)*t + s8)*t + s7)*t + s6) &
                            *t +  s5)*t +  s4)*t + s3)*t + s2)*t + s1)*t + s0
          erfc0 = rpinv*(z/x)
      endif

      end function erfc0

!*******************************************************************************
! exparg(l) =
!    if l == 0: the largest positive w for which exp(w) can be computed.
!    if l /= 0: the largest negative w for which the computed value of
!               exp(w) is nonzero.
!
! Note: only an approximate value for exparg(l) is needed.
!*******************************************************************************

      function exparg(l)

      use precision, only: wp

      implicit none
      integer, intent(in)  :: l
      real(wp)             :: exparg
      real(wp), parameter  :: one = 1.0_wp


      if (l == 0) then
          exparg = log(huge(one))
      else
          exparg = log(tiny(one))
      endif

      end function exparg

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

      end function erfc10
