!*******************************************************************************
! Complementary error function
! Taken from FNLIB/SLATEC
! july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
!*******************************************************************************

      function erfc12(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)  :: x
      real(wp)              :: erfc
      real(wp)              :: y,eta
      real(wp), parameter   :: sqrtpi = 1.77245385090551602729816748334115_wp
      real(wp), parameter   :: c1(21) = (/                                  &
                               -0.49046121234691808039984544033376e-01_wp,  &
                               -0.14226120510371364237824741899631e+00_wp,  &
                               +0.10035582187599795575754676712933e-01_wp,  &
                               -0.57687646997674847650827025509167e-03_wp,  &
                               +0.27419931252196061034422160791471e-04_wp,  &
                               -0.11043175507344507604135381295905e-05_wp,  &
                               +0.38488755420345036949961311498174e-07_wp,  &
                               -0.11808582533875466969631751801581e-08_wp,  &
                               +0.32334215826050909646402930953354e-10_wp,  &
                               -0.79910159470045487581607374708595e-12_wp,  &
                               +0.17990725113961455611967245486634e-13_wp,  &
                               -0.37186354878186926382316828209493e-15_wp,  &
                               +0.71035990037142529711689908394666e-17_wp,  &
                               -0.12612455119155225832495424853333e-18_wp,  &
                               +0.20916406941769294369170500266666e-20_wp,  &
                               -0.32539731029314072982364160000000e-22_wp,  &
                               +0.47668672097976748332373333333333e-24_wp,  &
                               -0.65980120782851343155199999999999e-26_wp,  &
                               +0.86550114699637626197333333333333e-28_wp,  &
                               -0.10788925177498064213333333333333e-29_wp,  &
                               +0.12811883993017002666666666666666e-31_wp  /)
      real(wp), parameter    :: c2(49) = (/                                 &
                                -0.6960134660230950112739150826197e-01_wp,  &
                                -0.4110133936262089348982212084666e-01_wp,  &
                                +0.3914495866689626881561143705244e-02_wp,  &
                                -0.4906395650548979161280935450774e-03_wp,  &
                                +0.7157479001377036380760894141825e-04_wp,  &
                                -0.1153071634131232833808232847912e-04_wp,  &
                                +0.1994670590201997635052314867709e-05_wp,  &
                                -0.3642666471599222873936118430711e-06_wp,  &
                                +0.6944372610005012589931277214633e-07_wp,  &
                                -0.1371220902104366019534605141210e-07_wp,  &
                                +0.2788389661007137131963860348087e-08_wp,  &
                                -0.5814164724331161551864791050316e-09_wp,  &
                                +0.1238920491752753181180168817950e-09_wp,  &
                                -0.2690639145306743432390424937889e-10_wp,  &
                                +0.5942614350847910982444709683840e-11_wp,  &
                                -0.1332386735758119579287754420570e-11_wp,  &
                                +0.3028046806177132017173697243304e-12_wp,  &
                                -0.6966648814941032588795867588954e-13_wp,  &
                                +0.1620854541053922969812893227628e-13_wp,  &
                                -0.3809934465250491999876913057729e-14_wp,  &
                                +0.9040487815978831149368971012975e-15_wp,  &
                                -0.2164006195089607347809812047003e-15_wp,  &
                                +0.5222102233995854984607980244172e-16_wp,  &
                                -0.1269729602364555336372415527780e-16_wp,  &
                                +0.3109145504276197583836227412951e-17_wp,  &
                                -0.7663762920320385524009566714811e-18_wp,  &
                                +0.1900819251362745202536929733290e-18_wp,  &
                                -0.4742207279069039545225655999965e-19_wp,  &
                                +0.1189649200076528382880683078451e-19_wp,  &
                                -0.3000035590325780256845271313066e-20_wp,  &
                                +0.7602993453043246173019385277098e-21_wp,  &
                                -0.1935909447606872881569811049130e-21_wp,  &
                                +0.4951399124773337881000042386773e-22_wp,  &
                                -0.1271807481336371879608621989888e-22_wp,  &
                                +0.3280049600469513043315841652053e-23_wp,  &
                                -0.8492320176822896568924792422399e-24_wp,  &
                                +0.2206917892807560223519879987199e-24_wp,  &
                                -0.5755617245696528498312819507199e-25_wp,  &
                                +0.1506191533639234250354144051199e-25_wp,  &
                                -0.3954502959018796953104285695999e-26_wp,  &
                                +0.1041529704151500979984645051733e-26_wp,  &
                                -0.2751487795278765079450178901333e-27_wp,  &
                                +0.7290058205497557408997703680000e-28_wp,  &
                                -0.1936939645915947804077501098666e-28_wp,  &
                                +0.5160357112051487298370054826666e-29_wp,  &
                                -0.1378419322193094099389644800000e-29_wp,  &
                                +0.3691326793107069042251093333333e-30_wp,  &
                                -0.9909389590624365420653226666666e-31_wp,  &
                                +0.2666491705195388413323946666666e-31_wp  /)
      real(wp), parameter    :: c3(59) = (/                                &
                                +0.715179310202924774503697709496e-01_wp,  &
                                -0.265324343376067157558893386681e-01_wp,  &
                                +0.171115397792085588332699194606e-02_wp,  &
                                -0.163751663458517884163746404749e-03_wp,  &
                                +0.198712935005520364995974806758e-04_wp,  &
                                -0.284371241276655508750175183152e-05_wp,  &
                                +0.460616130896313036969379968464e-06_wp,  &
                                -0.822775302587920842057766536366e-07_wp,  &
                                +0.159214187277090112989358340826e-07_wp,  &
                                -0.329507136225284321486631665072e-08_wp,  &
                                +0.722343976040055546581261153890e-09_wp,  &
                                -0.166485581339872959344695966886e-09_wp,  &
                                +0.401039258823766482077671768814e-10_wp,  &
                                -0.100481621442573113272170176283e-10_wp,  &
                                +0.260827591330033380859341009439e-11_wp,  &
                                -0.699111056040402486557697812476e-12_wp,  &
                                +0.192949233326170708624205749803e-12_wp,  &
                                -0.547013118875433106490125085271e-13_wp,  &
                                +0.158966330976269744839084032762e-13_wp,  &
                                -0.472689398019755483920369584290e-14_wp,  &
                                +0.143587337678498478672873997840e-14_wp,  &
                                -0.444951056181735839417250062829e-15_wp,  &
                                +0.140481088476823343737305537466e-15_wp,  &
                                -0.451381838776421089625963281623e-16_wp,  &
                                +0.147452154104513307787018713262e-16_wp,  &
                                -0.489262140694577615436841552532e-17_wp,  &
                                +0.164761214141064673895301522827e-17_wp,  &
                                -0.562681717632940809299928521323e-18_wp,  &
                                +0.194744338223207851429197867821e-18_wp,  &
                                -0.682630564294842072956664144723e-19_wp,  &
                                +0.242198888729864924018301125438e-19_wp,  &
                                -0.869341413350307042563800861857e-20_wp,  &
                                +0.315518034622808557122363401262e-20_wp,  &
                                -0.115737232404960874261239486742e-20_wp,  &
                                +0.428894716160565394623737097442e-21_wp,  &
                                -0.160503074205761685005737770964e-21_wp,  &
                                +0.606329875745380264495069923027e-22_wp,  &
                                -0.231140425169795849098840801367e-22_wp,  &
                                +0.888877854066188552554702955697e-23_wp,  &
                                -0.344726057665137652230718495566e-23_wp,  &
                                +0.134786546020696506827582774181e-23_wp,  &
                                -0.531179407112502173645873201807e-24_wp,  &
                                +0.210934105861978316828954734537e-24_wp,  &
                                -0.843836558792378911598133256738e-25_wp,  &
                                +0.339998252494520890627359576337e-25_wp,  &
                                -0.137945238807324209002238377110e-25_wp,  &
                                +0.563449031183325261513392634811e-26_wp,  &
                                -0.231649043447706544823427752700e-26_wp,  &
                                +0.958446284460181015263158381226e-27_wp,  &
                                -0.399072288033010972624224850193e-27_wp,  &
                                +0.167212922594447736017228709669e-27_wp,  &
                                -0.704599152276601385638803782587e-28_wp,  &
                                +0.297976840286420635412357989444e-28_wp,  &
                                -0.126252246646061929722422632994e-28_wp,  &
                                +0.539543870454248793985299653154e-29_wp,  &
                                -0.238099288253145918675346190062e-29_wp,  &
                                +0.109905283010276157359726683750e-29_wp,  &
                                -0.486771374164496572732518677435e-30_wp,  &
                                +0.152587726411035756763200828211e-30_wp  /)
      logical, save          :: first = .true.
      integer, save          :: n1 = 0
      integer, save          :: n3 = 0
      integer, save          :: n2 = 0
      real(wp), save         :: xmin = 0.0_wp
      real(wp), save         :: xmax = 0.0_wp
      real(wp), save         :: sqeps = 0.0_wp


      if (first) then
          first = .false.
          eta = 0.1*(r1mach(3))
          n1 = inits(c1,size(c1),eta)
          n3 = inits(c3,size(c3),eta)
          n2 = inits(c2,size(c2),eta)
          xmin = -sqrt(-log(sqrtpi*r1mach(3)))
          xmax = sqrt(-log(sqrtpi*r1mach(1)) )
          xmax = xmax - 0.5*log(xmax)/xmax - 0.01_wp
          sqeps = sqrt(2.0*r1mach(3))
      endif

      if (x <= xmin) then
          erfc = 2.0
      elseif (x > xmax) then
          erfc = 0.0
      else
          y = abs(x)
          if (y <= 1.0) then
              if (y < sqeps) then
                  erfc = 1.0 - 2.0*x/sqrtpi
              else
                  erfc = 1.0 - x*(1.0 + csevl(2.0*x*x-1.0,c1,n1))
              endif
          else
              y = y*y
              erfc = exp(-y)/abs(x)
              if (y <= 4.0) then
                  erfc = erfc * (0.5 + csevl((8.0/y-5.0)/3.0,c2,n2))
              else
                  erfc = erfc * (0.5 + csevl(8.0/y-1.0,c3,n3))
              endif
              if (x < 0.0) erfc = 2.0 - erfc
          endif
      endif

      contains

!*******************************************************************************
! Floating point machine constants (b=2)
! r1mach(1) = b**(emin-1), the smallest positive magnitude.
! r1mach(2) = b**emax*(1 - b**(-t)), the largest magnitude.
! r1mach(3) = b**(-t), the smallest relative spacing.
! r1mach(4) = b**(1-t), the largest relative spacing.
! r1mach(5) = log10(b)              
!*******************************************************************************

      function r1mach(i)

      use precision, only: wp

      implicit none
      integer, intent(in)       :: i
      real(wp)                  :: r1mach
      real(wp), save            :: rmach(5)
      logical, save             :: first = .true.
      real(wp), parameter       :: one = 1.0_wp


      if (first) then
          rmach(1) = tiny(one)
          rmach(2) = huge(one)
          rmach(3) = real(radix(one),wp)**(-digits(one))
          rmach(4) = epsilon(one)
          rmach(5) = log10(real(radix(one),wp))
          first = .false.
      endif

      select case(i)
      case(1:5)
          r1mach = rmach(i)
      case default
          stop 'R1MACH: i out of bounds'
      endselect

      end function r1mach

!*******************************************************************************
! april 1977 version.  w. fullerton, c3, los alamos scientific lab.
!
! initialize the orthogonal series so that inits is the number of terms
! needed to insure the error is no larger than eta.  ordinarily, eta
! will be chosen to be one-tenth machine precision.
!
! os     array of nos coefficients in an orthogonal series.
! nos    number of coefficients in os.
! eta    requested accuracy of series.
!*******************************************************************************

      function inits(os,nos,eta)

      use precision, only: wp

      implicit none
      integer, intent(in)        :: nos
      real(wp), intent(in)       :: os(nos),eta
      real(wp)                   :: err
      integer                    :: i,ii,inits


      if (nos < 1) stop 'INITS  nos < 1'

      err = 0.0_wp
      do ii=1,nos
          i = nos + 1 - ii
          err = err + abs(os(i))
          if (err > eta) exit
      enddo
      if (i == nos) stop 'INITS  eta may be too small'
      inits = i

      end function inits

!*******************************************************************************
! april 1977 version.  w. fullerton, c3, los alamos scientific lab.
!
! evaluate the n-term chebyshev series cs at x.  adapted from
! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).  also see fox
! and parker, chebyshev polys in numerical analysis, oxford press, p.56.
!
! x      value at which the series is to be evaluated.
! cs     array of n terms of a chebyshev series.  in eval-
!        uating cs, only half the first coef is summed.
! n      number of terms in array cs.
!*******************************************************************************

      function csevl(x,cs,n)

      use precision, only: wp

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


      if (n < 1) stop 'CSEVL n < 1'
      if (n > 1000) stop 'CSEVL n > 1000'
      if (x < -1.1_wp .or. x > 1.1_wp) stop 'CSEVL x outside (-1,+1)'

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

      csevl = (b0-b2)/2.0_wp

      end function csevl

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

      end function erfc12
