module eierf

  use const_and_precisions, only : wp_, zero, one
  implicit none
  real(wp_), parameter, private :: half=0.5_wp_, two=2.0_wp_, three=3.0_wp_, &
        four=4.0_wp_, six=6.0_wp_, twelve=12._wp_, sixten=16.0_wp_, &
        two4=24.0_wp_, fourty=40.0_wp_

contains

! ======================================================================
! nist guide to available math software.
! fullsource for module ei from package specfun.
! retrieved from netlib on fri mar 26 05:52:39 1999.
! ======================================================================
  subroutine calcei(arg,result,intt)
!----------------------------------------------------------------------
!
! this fortran 77 packet computes the exponential integrals ei(x),
!  e1(x), and  exp(-x)*ei(x)  for real arguments  x  where
!
!       integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
!  ei(x) =
!      -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
!
!  and where the first integral is a principal value integral.
!  the packet contains three function type subprograms: ei, eone,
!  and expei;  and one subroutine type subprogram: calcei.  the
!  calling statements for the primary entries are
!
!             y = ei(x),            where  x /= 0,
!
!             y = eone(x),          where  x > 0,
!  and
!             y = expei(x),         where  x /= 0,
!
!  and where the entry points correspond to the functions ei(x),
!  e1(x), and exp(-x)*ei(x), respectively.  the routine calcei
!  is intended for internal packet use only, all computations within
!  the packet being concentrated in this routine.  the function
!  subprograms invoke calcei with the fortran statement
!     call calcei(arg,result,intt)
!  where the parameter usage is as follows
!
! function                  parameters for calcei
!   call                 arg             result         intt
!
!  ei(x)              x /= 0          ei(x)            1
!  eone(x)            x > 0          -ei(-x)           2
!  expei(x)           x /= 0          exp(-x)*ei(x)    3
!----------------------------------------------------------------------
    implicit none
    integer, intent(in) :: intt
    real(wp_), intent(in) :: arg
    real(wp_), intent(out) :: result
    integer :: i
    real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq
    real(wp_), dimension(10) :: px,qx
!----------------------------------------------------------------------
!  mathematical constants
!   exp40 = exp(40)
!   x0 = zero of ei
!   x01/x11 + x02 = zero of ei to extra precision
!----------------------------------------------------------------------
    real(wp_), parameter :: p037=0.037_wp_, &
         exp40=2.3538526683701998541e17_wp_, x01=381.5_wp_, x11=1024.0_wp_, &
         x02=-5.1182968633365538008e-5_wp_, x0=3.7250741078136663466e-1_wp_
!----------------------------------------------------------------------
! machine-dependent constants
!----------------------------------------------------------------------
    real(wp_), parameter :: xinf=1.79e+308_wp_,xmax=716.351_wp_,xbig=701.84_wp_
!----------------------------------------------------------------------
! coefficients  for -1.0 <= x < 0.0
!----------------------------------------------------------------------
    real(wp_), dimension(7), parameter :: &
       a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, &
           1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, &
           1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, &
           5.0196785185439843791020_wp_/)
    real(wp_), dimension(6), parameter :: &
       b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, &
           8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, &
           1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/)
!----------------------------------------------------------------------
! coefficients for -4.0 <= x < -1.0
!----------------------------------------------------------------------
    real(wp_), dimension(9), parameter :: &
       c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, &
           7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, &
           1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, &
           1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, &
           1.737331760720576030932e-8_wp_/), &
       d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, &
           4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, &
           2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, &
           9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, &
           1.000000000000000000000e+0_wp_/)
!----------------------------------------------------------------------
! coefficients for x < -4.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, &
           1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, &
           1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, &
           1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, &
           5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/),&
       f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, &
           5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, &
           2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, &
           1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, &
           5.4199632588522559414924e+1_wp_,1.0_wp_/)
!----------------------------------------------------------------------
!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
!----------------------------------------------------------------------
    real(wp_), dimension(4), parameter :: &
       plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, &
             -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), &
       qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, &
             -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/)
!----------------------------------------------------------------------
! coefficients for  0.0 < x < 6.0,
!  ratio of chebyshev polynomials
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, &
           -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, &
           -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, &
            3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, &
            1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),&
       q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, &
            1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, &
            6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, &
            5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, &
            9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for 6.0 <= x < 12.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, &
           -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, &
            1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, &
            1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, &
            9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/)
    real(wp_), dimension(9), parameter :: &
       s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, &
            3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, &
            7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, &
            4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, &
            1.249884822712447891440_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for 12.0 <= x < 24.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, &
            -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, &
            -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, &
             2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, &
            -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/)
    real(wp_), dimension(9), parameter :: &
       q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, &
             5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, &
             4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, &
             1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, &
             1.001533852045342697818_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for  x >= 24.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, &
            -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, &
            -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, &
            -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, &
            -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/)
    real(wp_), dimension(9), parameter :: &
       q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, &
             1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, &
             7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, &
            -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, &
             1.99999999999048104167_wp_/)
!----------------------------------------------------------------------
    x = arg
    if (x == zero) then
      ei = -xinf
      if (intt == 2) ei = -ei
    else if ((x < zero) .or. (intt == 2)) then
!----------------------------------------------------------------------
! calculate ei for negative argument or for e1.
!----------------------------------------------------------------------
      y = abs(x)
      if (y <= one) then
        sump = a(7) * y + a(1)
        sumq = y + b(1)
        do i = 2, 6
          sump = sump * y + a(i)
          sumq = sumq * y + b(i)
        end do
        ei = log(y) - sump / sumq
        if (intt == 3) ei = ei * exp(y)
      else if (y <= four) then
        w = one / y
        sump = c(1)
        sumq = d(1)
        do i = 2, 9
          sump = sump * w + c(i)
          sumq = sumq * w + d(i)
        end do
        ei = - sump / sumq
        if (intt /= 3) ei = ei * exp(-y)
      else
        if ((y > xbig) .and. (intt < 3)) then
          ei = zero
        else
          w = one / y
          sump = e(1)
          sumq = f(1)
          do i = 2, 10
            sump = sump * w + e(i)
            sumq = sumq * w + f(i)
          end do
          ei = -w * (one - w * sump / sumq )
          if (intt /= 3) ei = ei * exp(-y)
        end if
      end if
      if (intt == 2) ei = -ei
    else if (x < six) then
!----------------------------------------------------------------------
!  to improve conditioning, rational approximations are expressed
!    in terms of chebyshev polynomials for 0 <= x < 6, and in
!    continued fraction form for larger x.
!----------------------------------------------------------------------
      t = x + x
      t = t / three - two
      px(1) = zero
      qx(1) = zero
      px(2) = p(1)
      qx(2) = q(1)
      do i = 2, 9
        px(i+1) = t * px(i) - px(i-1) + p(i)
        qx(i+1) = t * qx(i) - qx(i-1) + q(i)
      end do
      sump = half * t * px(10) - px(9) + p(10)
      sumq = half * t * qx(10) - qx(9) + q(10)
      frac = sump / sumq
      xmx0 = (x - x01/x11) - x02
      if (abs(xmx0) >= p037) then
        ei = log(x/x0) + xmx0 * frac
        if (intt == 3) ei = exp(-x) * ei
      else
!----------------------------------------------------------------------
! special approximation to  ln(x/x0)  for x close to x0
!----------------------------------------------------------------------
        y = xmx0 / (x + x0)
        ysq = y*y
        sump = plg(1)
        sumq = ysq + qlg(1)
        do i = 2, 4
          sump = sump*ysq + plg(i)
          sumq = sumq*ysq + qlg(i)
        end do
        ei = (sump / (sumq*(x+x0)) + frac) * xmx0
        if (intt == 3) ei = exp(-x) * ei
      end if
    else if (x < twelve) then
      frac = zero
      do i = 1, 9
        frac = s(i) / (r(i) + x + frac)
      end do
      ei = (r(10) + frac) / x
      if (intt /= 3) ei = ei * exp(x)
    else if (x <= two4) then
      frac = zero
      do i = 1, 9
        frac = q1(i) / (p1(i) + x + frac)
      end do
      ei = (p1(10) + frac) / x
      if (intt /= 3) ei = ei * exp(x)
    else
      if ((x >= xmax) .and. (intt < 3)) then
        ei = xinf
      else
        y = one / x
        frac = zero
        do i = 1, 9
          frac = q2(i) / (p2(i) + x + frac)
        end do
        frac = p2(10) + frac
        ei = y + y * y * frac
        if (intt /= 3) then
          if (x <= xmax-two4) then
            ei = ei * exp(x)
          else
!----------------------------------------------------------------------
! calculation reformulated to avoid premature overflow
!----------------------------------------------------------------------
            ei = (ei * exp(x-fourty)) * exp40
          end if
        end if
      end if
    end if
    result = ei
  end subroutine calcei

  function ei(x)
!--------------------------------------------------------------------
!
! this function program computes approximate values for the
!   exponential integral  ei(x), where  x  is real.
!
!  author: w. j. cody
!
!  latest modification: january 12, 1988
!
!--------------------------------------------------------------------
    implicit none
    integer :: intt
    real(wp_) :: ei
    real(wp_), intent(in) :: x
    real(wp_) :: result
!--------------------------------------------------------------------
    intt = 1
    call calcei(x,result,intt)
    ei = result
  end function ei

  function expei(x)
!--------------------------------------------------------------------
!
! this function program computes approximate values for the
!   function  exp(-x) * ei(x), where  ei(x)  is the exponential
!   integral, and  x  is real.
!
!  author: w. j. cody
!
!  latest modification: january 12, 1988
!
!--------------------------------------------------------------------
    implicit none
    integer :: intt
    real(wp_) :: expei
    real(wp_), intent(in) :: x
    real(wp_) :: result
!--------------------------------------------------------------------
    intt = 3
    call calcei(x,result,intt)
    expei = result
  end function expei

  function eone(x)
!--------------------------------------------------------------------
!
! this function program computes approximate values for the
!   exponential integral e1(x), where  x  is real.
!
!  author: w. j. cody
!
!  latest modification: january 12, 1988
!
!--------------------------------------------------------------------
    implicit none
    integer :: intt
    real(wp_) :: eone
    real(wp_), intent(in) :: x
    real(wp_) :: result
!--------------------------------------------------------------------
    intt = 2
    call calcei(x,result,intt)
    eone = result
  end function eone

! ======================================================================
! calcei3 = calcei for int=3
! ======================================================================
  subroutine calcei3(arg,result)
!----------------------------------------------------------------------
!
! this fortran 77 packet computes the exponential integrals ei(x),
!  e1(x), and  exp(-x)*ei(x)  for real arguments  x  where
!
!       integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
!  ei(x) =
!      -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
!
!  and where the first integral is a principal value integral.
!  the packet contains three function type subprograms: ei, eone,
!  and expei;  and one subroutine type subprogram: calcei.  the
!  calling statements for the primary entries are
!
!             y = ei(x),            where  x /= 0,
!
!             y = eone(x),          where  x > 0,
!  and
!             y = expei(x),         where  x /= 0,
!
!  and where the entry points correspond to the functions ei(x),
!  e1(x), and exp(-x)*ei(x), respectively.  the routine calcei
!  is intended for internal packet use only, all computations within
!  the packet being concentrated in this routine.  the function
!  subprograms invoke calcei with the fortran statement
!     call calcei(arg,result,int)
!  where the parameter usage is as follows
!
! function                  parameters for calcei
!   call                 arg             result         int
!
!  ei(x)              x /= 0          ei(x)            1
!  eone(x)            x > 0         -ei(-x)           2
!  expei(x)           x /= 0          exp(-x)*ei(x)    3
!----------------------------------------------------------------------
    implicit none
    real(wp_), intent(in) :: arg
    real(wp_), intent(out) :: result
    integer :: i
    real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq
    real(wp_), dimension(10) :: px,qx
!----------------------------------------------------------------------
!  mathematical constants
!   exp40 = exp(40)
!   x0 = zero of ei
!   x01/x11 + x02 = zero of ei to extra precision
!----------------------------------------------------------------------
    real(wp_), parameter :: p037=0.037_wp_, &
         x01=381.5_wp_, x11=1024.0_wp_, x02=-5.1182968633365538008e-5_wp_, &
         x0=3.7250741078136663466e-1_wp_
!----------------------------------------------------------------------
! machine-dependent constants
!----------------------------------------------------------------------
    real(wp_), parameter :: xinf=1.79e+308_wp_
!----------------------------------------------------------------------
! coefficients  for -1.0 <= x < 0.0
!----------------------------------------------------------------------
    real(wp_), dimension(7), parameter :: &
       a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, &
           1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, &
           1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, &
           5.0196785185439843791020_wp_/)
    real(wp_), dimension(6), parameter :: &
       b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, &
           8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, &
           1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/)
!----------------------------------------------------------------------
! coefficients for -4.0 <= x < -1.0
!----------------------------------------------------------------------
    real(wp_), dimension(9), parameter :: &
       c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, &
           7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, &
           1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, &
           1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, &
           1.737331760720576030932e-8_wp_/), &
       d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, &
           4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, &
           2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, &
           9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, &
           1.000000000000000000000e+0_wp_/)
!----------------------------------------------------------------------
! coefficients for x < -4.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, &
           1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, &
           1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, &
           1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, &
           5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/), &
       f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, &
           5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, &
           2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, &
           1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, &
           5.4199632588522559414924e+1_wp_,1.0_wp_/)
!----------------------------------------------------------------------
!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
!----------------------------------------------------------------------
    real(wp_), dimension(4), parameter :: &
       plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, &
             -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), &
       qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, &
             -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/)
!----------------------------------------------------------------------
! coefficients for  0.0 < x < 6.0,
!  ratio of chebyshev polynomials
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, &
           -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, &
           -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, &
            3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, &
            1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),&
       q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, &
            1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, &
            6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, &
            5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, &
            9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for 6.0 <= x < 12.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, &
           -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, &
            1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, &
            1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, &
            9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/)
    real(wp_), dimension(9), parameter :: &
       s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, &
            3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, &
            7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, &
            4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, &
            1.249884822712447891440_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for 12.0 <= x < 24.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, &
            -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, &
            -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, &
             2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, &
            -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/)
    real(wp_), dimension(9), parameter :: &
       q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, &
             5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, &
             4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, &
             1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, &
             1.001533852045342697818_wp_/)
!----------------------------------------------------------------------
! j-fraction coefficients for  x >= 24.0
!----------------------------------------------------------------------
    real(wp_), dimension(10), parameter :: &
       p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, &
            -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, &
            -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, &
            -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, &
            -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/)
    real(wp_), dimension(9), parameter :: &
       q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, &
             1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, &
             7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, &
            -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, &
             1.99999999999048104167_wp_/)
!----------------------------------------------------------------------
    x = arg
    if (x == zero) then
      ei = -xinf
    else if ((x < zero)) then
!----------------------------------------------------------------------
! calculate ei for negative argument or for e1.
!----------------------------------------------------------------------
      y = abs(x)
      if (y <= one) then
        sump = a(7) * y + a(1)
        sumq = y + b(1)
        do i = 2, 6
          sump = sump * y + a(i)
          sumq = sumq * y + b(i)
        end do
        ei = (log(y) - sump / sumq ) * exp(y)
      else if (y <= four) then
        w = one / y
        sump = c(1)
        sumq = d(1)
        do i = 2, 9
          sump = sump * w + c(i)
          sumq = sumq * w + d(i)
        end do
        ei = - sump / sumq
      else
        w = one / y
        sump = e(1)
        sumq = f(1)
        do i = 2, 10
          sump = sump * w + e(i)
          sumq = sumq * w + f(i)
        end do
        ei = -w * (one - w * sump / sumq )
      end if
    else if (x < six) then
!----------------------------------------------------------------------
!  to improve conditioning, rational approximations are expressed
!    in terms of chebyshev polynomials for 0 <= x < 6, and in
!    continued fraction form for larger x.
!----------------------------------------------------------------------
      t = x + x
      t = t / three - two
      px(1) = zero
      qx(1) = zero
      px(2) = p(1)
      qx(2) = q(1)
      do i = 2, 9
        px(i+1) = t * px(i) - px(i-1) + p(i)
        qx(i+1) = t * qx(i) - qx(i-1) + q(i)
      end do
      sump = half * t * px(10) - px(9) + p(10)
      sumq = half * t * qx(10) - qx(9) + q(10)
      frac = sump / sumq
      xmx0 = (x - x01/x11) - x02
      if (abs(xmx0) >= p037) then
        ei = exp(-x) * ( log(x/x0) + xmx0 * frac )
      else
!----------------------------------------------------------------------
! special approximation to  ln(x/x0)  for x close to x0
!----------------------------------------------------------------------
        y = xmx0 / (x + x0)
        ysq = y*y
        sump = plg(1)
        sumq = ysq + qlg(1)
        do i = 2, 4
          sump = sump*ysq + plg(i)
          sumq = sumq*ysq + qlg(i)
        end do
        ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0
      end if
    else if (x < twelve) then
      frac = zero
      do i = 1, 9
        frac = s(i) / (r(i) + x + frac)
      end do
      ei = (r(10) + frac) / x
    else if (x <= two4) then
      frac = zero
      do i = 1, 9
        frac = q1(i) / (p1(i) + x + frac)
      end do
      ei = (p1(10) + frac) / x
    else
      y = one / x
      frac = zero
      do i = 1, 9
        frac = q2(i) / (p2(i) + x + frac)
      end do
      frac = p2(10) + frac
      ei = y + y * y * frac
    end if
    result = ei
  end subroutine calcei3

!  subroutine calerf(arg,result,jintt)
!!------------------------------------------------------------------
!!
!! this packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x)
!!   for a real argument  x.  it contains three function type
!!   subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx),
!!   and one subroutine type subprogram, calerf.  the calling
!!   statements for the primary entries are:
!!
!!               y=erf(x)     (or   y=derf(x)),
!!
!!               y=erfc(x)    (or   y=derfc(x)),
!!   and
!!               y=erfcx(x)   (or   y=derfcx(x)).
!!
!!   the routine  calerf  is intended for internal packet use only,
!!   all computations within the packet being concentrated in this
!!   routine.  the function subprograms invoke  calerf  with the
!!   statement
!!
!!      call calerf(arg,result,jintt)
!!
!!   where the parameter usage is as follows
!!
!!  function                     parameters for calerf
!!   call              arg                  result          jintt
!!
!! erf(arg)      any real argument         erf(arg)          0
!! erfc(arg)     abs(arg) < xbig           erfc(arg)         1
!! erfcx(arg)    xneg < arg < xmax         erfcx(arg)        2
!!
!!*******************************************************************
!!*******************************************************************
!!
!! Explanation of machine-dependent constants
!!
!!   XMIN   = the smallest positive floating-point number.
!!   XINF   = the largest positive finite floating-point number.
!!   XNEG   = the largest negative argument acceptable to ERFCX;
!!            the negative of the solution to the equation
!!            2*exp(x*x) = XINF.
!!   XSMALL = argument below which erf(x) may be represented by
!!            2*x/sqrt(pi)  and above which  x*x  will not underflow.
!!            A conservative value is the largest machine number X
!!            such that   1.0 + X = 1.0   to machine precision.
!!   XBIG   = largest argument acceptable to ERFC;  solution to
!!            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where
!!            W(x) = exp(-x*x)/[x*sqrt(pi)].
!!   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to
!!            machine precision.  A conservative value is
!!            1/[2*sqrt(XSMALL)]
!!   XMAX   = largest acceptable argument to ERFCX; the minimum
!!            of XINF and 1/[sqrt(pi)*XMIN].
!!
!!*******************************************************************
!!*******************************************************************
!!
!! error returns
!!
!!  the program returns  erfc = 0      for  arg >= xbig;
!!
!!                   erfcx = xinf  for  arg < xneg;
!!  and
!!                   erfcx = 0     for  arg >= xmax.
!!
!!
!! intrinsic functions required are:
!!
!! abs, aint, exp
!!
!!
!!  author: w. j. cody
!!      mathematics and computer science division
!!      argonne national laboratory
!!      argonne, il 60439
!!
!!  latest modification: march 19, 1990
!!
!!------------------------------------------------------------------
!    implicit none
!    real(wp_), intent(in) :: arg
!    real(wp_), intent(out) :: result
!    integer, intent(in) :: jintt
!    integer :: i
!    real(wp_) :: del,x,xden,xnum,y,ysq
!!------------------------------------------------------------------
!!  mathematical constants
!!------------------------------------------------------------------
!    real(wp_), parameter :: sqrpi=5.6418958354775628695e-1_wp_, &
!       thresh=0.46875_wp_
!!------------------------------------------------------------------
!!  machine-dependent constants
!!------------------------------------------------------------------
!    real(wp_), parameter :: xinf=1.79e308_wp_,    &  ! ~huge
!                            xneg=-26.628_wp_,     &  ! ?
!                            xsmall=1.11e-16_wp_,  &  ! ~epsilon/2
!                            xbig=26.543_wp_,      &  ! ?
!                            xhuge=6.71e7_wp_,     &  ! ~1/sqrt(epsilon)
!                            xmax=2.53e307_wp_        ! ?
!!------------------------------------------------------------------
!!  coefficients for approximation to  erf  in first interval
!!------------------------------------------------------------------
!    real(wp_), dimension(5), parameter :: &
!       a=(/3.16112374387056560_wp_,1.13864154151050156e02_wp_, &
!           3.77485237685302021e02_wp_,3.20937758913846947e03_wp_, &
!           1.85777706184603153e-1_wp_/)
!    real(wp_), dimension(4), parameter :: &
!       b=(/2.36012909523441209e01_wp_,2.44024637934444173e02_wp_, &
!           1.28261652607737228e03_wp_,2.84423683343917062e03_wp_/)
!!------------------------------------------------------------------
!!  coefficients for approximation to  erfc  in second interval
!!------------------------------------------------------------------
!    real(wp_), dimension(9), parameter :: &
!       c=(/5.64188496988670089e-1_wp_,8.88314979438837594_wp_, &
!           6.61191906371416295e01_wp_,2.98635138197400131e02_wp_, &
!           8.81952221241769090e02_wp_,1.71204761263407058e03_wp_, &
!           2.05107837782607147e03_wp_,1.23033935479799725e03_wp_, &
!           2.15311535474403846e-8_wp_/)
!    real(wp_), dimension(8), parameter :: &
!       d=(/1.57449261107098347e01_wp_,1.17693950891312499e02_wp_, &
!           5.37181101862009858e02_wp_,1.62138957456669019e03_wp_, &
!           3.29079923573345963e03_wp_,4.36261909014324716e03_wp_, &
!           3.43936767414372164e03_wp_,1.23033935480374942e03_wp_/)
!!------------------------------------------------------------------
!!  coefficients for approximation to  erfc  in third interval
!!------------------------------------------------------------------
!    real(wp_), dimension(6), parameter :: &
!       p=(/3.05326634961232344e-1_wp_,3.60344899949804439e-1_wp_, &
!           1.25781726111229246e-1_wp_,1.60837851487422766e-2_wp_, &
!           6.58749161529837803e-4_wp_,1.63153871373020978e-2_wp_/)
!    real(wp_), dimension(5), parameter :: &
!       q=(/2.56852019228982242_wp_,1.87295284992346047_wp_, &
!           5.27905102951428412e-1_wp_,6.05183413124413191e-2_wp_, &
!           2.33520497626869185e-3_wp_/)
!!------------------------------------------------------------------
!    x = arg
!    y = abs(x)
!    if (y <= thresh) then
!!------------------------------------------------------------------
!!  evaluate  erf  for  |x| <= 0.46875
!!------------------------------------------------------------------
!      ysq = zero
!      if (y > xsmall) ysq = y * y
!      xnum = a(5)*ysq
!      xden = ysq
!      do i = 1, 3
!        xnum = (xnum + a(i)) * ysq
!        xden = (xden + b(i)) * ysq
!      end do
!      result = x * (xnum + a(4)) / (xden + b(4))
!      if (jintt /= 0) result = one - result
!      if (jintt == 2) result = exp(ysq) * result
!      return
!!------------------------------------------------------------------
!!  evaluate  erfc  for 0.46875 <= |x| <= 4.0
!!------------------------------------------------------------------
!    else if (y <= four) then
!      xnum = c(9)*y
!      xden = y
!      do i = 1, 7
!        xnum = (xnum + c(i)) * y
!        xden = (xden + d(i)) * y
!      end do
!      result = (xnum + c(8)) / (xden + d(8))
!      if (jintt /= 2) then
!        ysq = aint(y*sixten)/sixten
!        del = (y-ysq)*(y+ysq)
!        result = exp(-ysq*ysq) * exp(-del) * result
!      end if
!!------------------------------------------------------------------
!!  evaluate  erfc  for |x| > 4.0
!!------------------------------------------------------------------
!    else if (y < xbig .or. (y < xmax .and. jintt == 2)) then
!      ysq = one / (y * y)
!      xnum = p(6)*ysq
!      xden = ysq
!      do i = 1, 4
!        xnum = (xnum + p(i)) * ysq
!        xden = (xden + q(i)) * ysq
!      end do
!      result = ysq *(xnum + p(5)) / (xden + q(5))
!      result = (sqrpi -  result) / y
!      if (jintt /= 2) then
!        ysq = aint(y*sixten)/sixten
!        del = (y-ysq)*(y+ysq)
!        result = exp(-ysq*ysq) * exp(-del) * result
!      end if
!    else if (y >= xhuge) then
!      result = sqrpi / y
!    else
!      result = zero
!    end if
!!------------------------------------------------------------------
!!  fix up for negative argument, erf, etc.
!!------------------------------------------------------------------
!    if (jintt == 0) then
!      result = (half - result) + half
!      if (x < zero) result = -result
!    else if (jintt == 1) then
!      if (x < zero) result = two - result
!    else
!      if (x < zero) then
!        if (x < xneg) then
!          result = xinf
!        else
!          ysq = aint(x*sixten)/sixten
!          del = (x-ysq)*(x+ysq)
!          y = exp(ysq*ysq) * exp(del)
!          result = (y+y) - result
!        end if
!      end if
!    end if
!  end subroutine calerf
!
!  function derf(x)
!!--------------------------------------------------------------------
!!
!! this subprogram computes approximate values for erf(x).
!!   (see comments heading calerf).
!!
!!   author/date: w. j. cody, january 8, 1985
!!
!!--------------------------------------------------------------------
!    implicit none
!    real(wp_) :: derf
!    real(wp_), intent(in) :: x
!    integer :: jintt
!    real(wp_) :: result
!!------------------------------------------------------------------
!    jintt = 0
!    call calerf(x,result,jintt)
!    derf = result
!  end function derf
!
!  function derfc(x)
!!--------------------------------------------------------------------
!!
!! this subprogram computes approximate values for erfc(x).
!!   (see comments heading calerf).
!!
!!   author/date: w. j. cody, january 8, 1985
!!
!!--------------------------------------------------------------------
!    implicit none
!    real(wp_) :: derfc
!    real(wp_), intent(in) :: x
!    integer :: jintt
!    real(wp_) :: result
!!------------------------------------------------------------------
!    jintt = 1
!    call calerf(x,result,jintt)
!    derfc = result
!  end function derfc
!
!  function derfcx(x)
!!------------------------------------------------------------------
!!
!! this subprogram computes approximate values for exp(x*x) * erfc(x).
!!   (see comments heading calerf).
!!
!!   author/date: w. j. cody, march 30, 1987
!!
!!------------------------------------------------------------------
!    implicit none
!    real(wp_) :: derfcx
!    real(wp_), intent(in) :: x
!    integer :: jintt
!    real(wp_) :: result
!!------------------------------------------------------------------
!    jintt = 2
!    call calerf(x,result,jintt)
!    derfcx = result
!  end function derfcx

end module eierf