anchor.diff:      last = 0
       result = 0.0_real64
       abserr = 0.0_real64
  -    alist(1) = a
  -    blist(1) = b
  +    alist(1) = 0.0_real64
  +    blist(1) = 0.1e1_real64
       rlist(1) = 0.0_real64
       elist(1) = 0.0_real64
       iord(1) = 0
   
  -    if(epsabs.le.0.0_real64.and. &
  -         epsrel.lt. max ( 0.5e2_real64*epmach,0.5e-28_real64)) then
  +    if(epsabs.le.0.0_real64.and.epsrel.lt. max ( 0.5e2_real64*epmach,0.5D-28)) then
          ier = 6
  +    end if
  +
  +    if(ier.eq.6) then
          return
       end if
       !
       !  first approximation to the integral
       !
  -    keyf = key
  -    if(key.le.0) keyf = 1
  -    if(key.ge.7) keyf = 6
  -    neval = 0
  -    if(keyf.eq.1) call dqk15(f,a,b,result,abserr,defabs,resabs)
  -    if(keyf.eq.2) call dqk21(f,a,b,result,abserr,defabs,resabs)
  -    if(keyf.eq.3) call dqk31(f,a,b,result,abserr,defabs,resabs)
  -    if(keyf.eq.4) call dqk41(f,a,b,result,abserr,defabs,resabs)
  -    if(keyf.eq.5) call dqk51(f,a,b,result,abserr,defabs,resabs)
  -    if(keyf.eq.6) call dqk61(f,a,b,result,abserr,defabs,resabs)
  +    !  determine the interval to be mapped onto (0,1).
  +    !  if inf = 2 the integral is computed as i = i1+i2, where
  +    !  i1 = integral of f over (-infinity,0),
  +    !  i2 = integral of f over (0,+infinity).
  +    !
  +    boun = bound
  +    if(inf.eq.2) boun = 0.0_real64
  +    call dqk15i(f,boun,inf,0.0_real64,0.1e1_real64,result,abserr, &
  +         defabs,resabs)
  +    !
  +    !  test on accuracy
  +    !
       last = 1
       rlist(1) = result
       elist(1) = abserr
       iord(1) = 1
  -    !
  -    !  test on accuracy.
  -    !
  -    errbnd =  max ( epsabs, epsrel* abs ( result ) )
  -
  -    if(abserr.le.0.5e2_real64 * epmach * defabs .and. &
  -         abserr.gt.errbnd) then
  -       ier = 2
  -    end if
  -
  -    if(limit.eq.1) then
  -       ier = 1
  -    end if
  -
  -    if ( ier .ne. 0 .or. &
  -         (abserr .le. errbnd .and. abserr .ne. resabs ) .or. &
  -         abserr .eq. 0.0_real64 ) then
  +    dres =  abs ( result)
  +    errbnd =  max ( epsabs,epsrel*dres)
  +    if(abserr.le.1.0e2_real64*epmach*defabs.and.abserr.gt.errbnd) ier = 2
  +    if(limit.eq.1) ier = 1
  +    if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or. &
  +         abserr.eq.0.0_real64) go to 130
  +    !
  +    !  initialization
  +    !
  +    uflow = tiny ( uflow )
  +    oflow = huge ( oflow )
  +    rlist2(1) = result
  +    errmax = abserr
  +    maxerr = 1
  +    area = result
  +    errsum = abserr
  +    abserr = oflow
  +    nrmax = 1
  +    nres = 0
  +    ktmin = 0
  +    numrl2 = 2
  +    extrap = .false.
  +    noext = .false.
  +    ierro = 0
  +    iroff1 = 0
  +    iroff2 = 0
  +    iroff3 = 0
  +    ksgn = -1
  +    if(dres.ge.(0.1e1_real64-0.5e2_real64*epmach)*defabs) ksgn = 1
  +    !
  +    !  main do-loop
  +    !
  +    do 90 last = 2,limit
  +       !
  +       !  bisect the subinterval with nrmax-th largest error estimate.
  +       !
  +       a1 = alist(maxerr)
  +       b1 = 0.5_real64*(alist(maxerr)+blist(maxerr))
  +       a2 = b1
  +       b2 = blist(maxerr)
  +       erlast = errmax
  +       call dqk15i(f,boun,inf,a1,b1,area1,error1,resabs,defab1)
  +       call dqk15i(f,boun,inf,a2,b2,area2,error2,resabs,defab2)
  +       !
  +       !  improve previous approximations to integral
  +       !  and error and test for accuracy.
  +       !
  +       area12 = area1+area2
  +       erro12 = error1+error2
  +       errsum = errsum+erro12-errmax
  +       area = area+area12-rlist(maxerr)
  +       if(defab1.eq.error1.or.defab2.eq.error2)go to 15
  +       if( abs ( rlist(maxerr)-area12).gt.0.1e-4_real64* abs ( area12) &
  +            .or.erro12.lt.0.99_real64*errmax) go to 10
  +       if(extrap) iroff2 = iroff2+1
  +       if(.not.extrap) iroff1 = iroff1+1
  +10     if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
  +15     rlist(maxerr) = area1
  +       rlist(last) = area2
  +       errbnd =  max ( epsabs,epsrel* abs ( area))
  +       !
  +       !  test for roundoff error and eventually set error flag.
  +       !
  +       if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
  +       if(iroff2.ge.5) ierro = 3
  +       !
  +       !  set error flag in the case that the number of
  +       !  subintervals equals limit.
  +       !
  +       if(last.eq.limit) ier = 1
  +       !
  +       !  set error flag in the case of bad integrand behaviour
  +       !  at some points of the integration range.
  +       !
  +       if( max (  abs ( a1), abs ( b2)).le.(0.1e1_real64+0.1e3_real64*epmach)* &
  +            ( abs ( a2)+0.1e4_real64*uflow)) then
  +          ier = 4
  +       end if
  +       !
  +       !  append the newly-created intervals to the list.
  +       !
  +       if(error2.gt.error1) go to 20
  +       alist(last) = a2
  +       blist(maxerr) = b1
  +       blist(last) = b2
  +       elist(maxerr) = error1
  +       elist(last) = error2
  +       go to 30
  +20     continue
  +
  +       alist(maxerr) = a2
  +       alist(last) = a1
  +       blist(last) = b1
  +       rlist(maxerr) = area2
  +       rlist(last) = area1
  +       elist(maxerr) = error2
  +       elist(last) = error1
  +       !
  +       !  call dqpsrt to maintain the descending ordering
  +       !  in the list of error estimates and select the subinterval
  +       !  with nrmax-th largest error estimate (to be bisected next).
  +       !
  +30     call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
  +       if(errsum.le.errbnd) go to 115
  +       if(ier.ne.0) go to 100
  +       if(last.eq.2) go to 80
  +       if(noext) go to 90
  +       erlarg = erlarg-erlast
  +       if( abs ( b1-a1).gt.small) erlarg = erlarg+erro12
  +       if(extrap) go to 40
  +       !
  +       !  test whether the interval to be bisected next is the
  +       !  smallest interval.
  +       !
  +       if( abs ( blist(maxerr)-alist(maxerr)).gt.small) go to 90
  +       extrap = .true.
  +       nrmax = 2
  +40     if(ierro.eq.3.or.erlarg.le.ertest) go to 60
  +       !
  +       !  the smallest interval has the largest error.
  +       !  before bisecting decrease the sum of the errors over the
  +       !  larger intervals (erlarg) and perform extrapolation.
  +       !
  +       id = nrmax
  +       jupbnd = last
  +       if(last.gt.(2+limit/2)) jupbnd = limit+3-last
  +
  +       do k = id,jupbnd
  +          maxerr = iord(nrmax)
  +          errmax = elist(maxerr)
  +          if( abs ( blist(maxerr)-alist(maxerr)).gt.small) go to 90
  +          nrmax = nrmax+1
  +       end do
  +       !
  +       !  perform extrapolation.
  +       !
  +60     numrl2 = numrl2+1
  +       rlist2(numrl2) = area
  +       call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
  +       ktmin = ktmin+1
  +       if(ktmin.gt.5.and.abserr.lt.0.1e-2_real64*errsum) ier = 5
  +       if(abseps.ge.abserr) go to 70
  +       ktmin = 0
  +       abserr = abseps
  +       result = reseps
  +       correc = erlarg
  +       ertest =  max ( epsabs,epsrel* abs ( reseps))
  +       if(abserr.le.ertest) go to 100
  +       !
  +       !  prepare bisection of the smallest interval.
  +       !
  +70     if(numrl2.eq.1) noext = .true.
  +       if(ier.eq.5) go to 100
  +       maxerr = iord(1)
  +       errmax = elist(maxerr)
  +       nrmax = 1
  +       extrap = .false.
  +       small = small*0.5_real64
  +       erlarg = errsum
  +       go to 90
  +80     small = 0.375_real64
  +       erlarg = errsum
  +       ertest = errbnd
  +       rlist2(2) = area
  +90     continue
  +       !
  +       !  set final result and error estimate.
  +       !
  +100    if(abserr.eq.oflow) go to 115
  +       if((ier+ierro).eq.0) go to 110
  +       if(ierro.eq.3) abserr = abserr+correc
  +       if(ier.eq.0) ier = 3
  +       if(result.ne.0.0_real64.and.area.ne.0.0_real64)go to 105
  +       if(abserr.gt.errsum)go to 115
  +       if(area.eq.0.0_real64) go to 130
  +       go to 110
  +105    if(abserr/ abs ( result).gt.errsum/ abs ( area))go to 115
  +       !
  +       !  test on divergence
  +       !
  +110    continue
  +
  +       if ( ksgn .eq. (-1) .and. &
  +            max ( abs ( result), abs ( area)) .le. defabs*0.1e-1_real64 ) then
  +          go to 130
  +       end if
  +
  +       if ( 0.1e-1_real64 .gt. (result/area) .or. &
  +            (result/area) .gt. 0.1e3_real64 .or. &
  +            errsum .gt. abs ( area) ) then
  +          ier = 6
  +       end if
  +
  +       go to 130
  +       !
  +       !  compute global integral sum.
  +       !
  +115    result = 0.0_real64
  +       do k = 1,last
  +          result = result+rlist(k)
  +       end do
  +       abserr = errsum
  +130    continue
  +
  +       neval = 30*last-15
  +       if(inf.eq.2) neval = 2*neval
  +       if(ier.gt.2) ier=ier-1
  +
  +       return
  +     end subroutine dqagie
  +     subroutine dqagi ( f, bound, inf, epsabs, epsrel, result, abserr, neval, &
  +          ier,limit,lenw,last,iwork,work)
  +
  +       !*****************************************************************************80
  +       !
  +       !! DQAGI estimates an integral over a semi-infinite or infinite interval.
  +       !
  +       !  Modified:
  +       !
  +       !    11 September 2015
  +       !
  +       !  Author:
  +       !
  +       !    Robert Piessens, Elise de Doncker
  +       !
  +       !***purpose  the routine calculates an approximation result to a given
  +       !      integral   i = integral of f over (bound,+infinity)
  +       !      or i = integral of f over (-infinity,bound)
  +       !      or i = integral of f over (-infinity,+infinity)
  +       !      hopefully satisfying following claim for accuracy
  +       !      abs(i-result).le.max(epsabs,epsrel*abs(i)).
  +       !
  +       !  Parameters:
  +       !
  +       !   on entry
  +       !      f      - real ( kind = 8 )
  +       !               function subprogram defining the integrand
  +       !               function f(x). the actual name for f needs to be
  +       !               declared e x t e r n a l in the driver program.
  +       !
  +       !      bound  - real ( kind = 8 )
  +       !               finite bound of integration range
  +       !               (has no meaning if interval is doubly-infinite)
  +       !
  +       !      inf    - integer ( kind = 4 )
  +       !               indicating the kind of integration range involved
  +       !               inf = 1 corresponds to  (bound,+infinity),
  +       !               inf = -1            to  (-infinity,bound),
  +       !               inf = 2             to (-infinity,+infinity).
  +       !
  +       !      epsabs - real ( kind = 8 )
  +       !               absolute accuracy requested
  +       !      epsrel - real ( kind = 8 )
  +       !               relative accuracy requested
  +       !               if  epsabs.le.0
  +       !               and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
  +       !               the routine will end with ier = 6.
  +       !
  +       !
  +       !   on return
  +       !      result - real ( kind = 8 )
  +       !               approximation to the integral
  +       !
  +       !      abserr - real ( kind = 8 )
  +       !               estimate of the modulus of the absolute error,
  +       !               which should equal or exceed abs(i-result)
  +       !
  +       !      neval  - integer ( kind = 4 )
  +       !               number of integrand evaluations
  +       !
  +       !      ier    - integer ( kind = 4 )
  +       !               ier = 0 normal and reliable termination of the
  +       !                       routine. it is assumed that the requested
  +       !                       accuracy has been achieved.
  +       !             - ier.gt.0 abnormal termination of the routine. the
  +       !                       estimates for result and error are less
  +       !                       reliable. it is assumed that the requested
  +       !                       accuracy has not been achieved.
  +       !      error messages
  +       !               ier = 1 maximum number of subdivisions allowed
  +       !                       has been achieved. one can allow more
  +       !                       subdivisions by increasing the value of
  +       !                       limit (and taking the according dimension
  +       !                       adjustments into account). however, if
  +       !                       this yields no improvement it is advised
  +       !                       to analyze the integrand in order to
  +       !                       determine the integration difficulties. if
  +       !                       the position of a local difficulty can be
  +       !                       determined (e.g. singularity,
  +       !                       discontinuity within the interval) one
  +       !                       will probably gain from splitting up the
  +       !                       interval at this point and calling the
  +       !                       integrator on the subranges. if possible,
  +       !                       an appropriate special-purpose integrator
  +       !                       should be used, which is designed for
  +       !                       handling the type of difficulty involved.
  +       !                   = 2 the occurrence of roundoff error is
  +       !                       detected, which prevents the requested
  +       !                       tolerance from being achieved.
  +       !                       the error may be under-estimated.
  +       !                   = 3 extremely bad integrand behaviour occurs
  +       !                       at some points of the integration
  +       !                       interval.
  +       !                   = 4 the algorithm does not converge.
  +       !                       roundoff error is detected in the
  +       !                       extrapolation table.
  +       !                       it is assumed that the requested tolerance
  +       !                       cannot be achieved, and that the returned
  +       !                       result is the best which can be obtained.
  +       !                   = 5 the integral is probably divergent, or
  +       !                       slowly convergent. it must be noted that
  +       !                       divergence can occur with any other value
  +       !                       of ier.
  +       !                   = 6 the input is invalid, because
  +       !                       (epsabs.le.0 and
  +       !                        epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
  +       !                        or limit.lt.1 or leniw.lt.limit*4.
  +       !                       result, abserr, neval, last are set to
  +       !                       zero. exept when limit or leniw is
  +       !                       invalid, iwork(1), work(limit*2+1) and
  +       !                       work(limit*3+1) are set to zero, work(1)
  +       !                       is set to a and work(limit+1) to b.
  +       !
  +       !   dimensioning parameters
  +       !      limit - integer ( kind = 4 )
  +       !              dimensioning parameter for iwork
  +       !              limit determines the maximum number of subintervals
  +       !              in the partition of the given integration interval
  +       !              (a,b), limit.ge.1.
  +       !              if limit.lt.1, the routine will end with ier = 6.
  +       !
  +       !      lenw  - integer ( kind = 4 )
  +       !              dimensioning parameter for work
  +       !              lenw must be at least limit*4.
  +       !              if lenw.lt.limit*4, the routine will end
  +       !              with ier = 6.
  +       !
  +       !      last  - integer ( kind = 4 )
  +       !              on return, last equals the number of subintervals
  +       !              produced in the subdivision process, which
  +       !              determines the number of significant elements
  +       !              actually in the work arrays.
  +       !
  +       !   work arrays
  +       !      iwork - integer ( kind = 4 )
  +       !              vector of dimension at least limit, the first
  +       !              k elements of which contain pointers
  +       !              to the error estimates over the subintervals,
  +       !              such that work(limit*3+iwork(1)),... ,
  +       !              work(limit*3+iwork(k)) form a decreasing
  +       !              sequence, with k = last if last.le.(limit/2+2), and
  +       !              k = limit+1-last otherwise
  +       !
  +       !      work  - real ( kind = 8 )
  +       !              vector of dimension at least lenw
  +       !              on return
  +       !              work(1), ..., work(last) contain the left
  +       !               end points of the subintervals in the
  +       !               partition of (a,b),
  +       !              work(limit+1), ..., work(limit+last) contain
  +       !               the right end points,
  +       !              work(limit*2+1), ...,work(limit*2+last) contain the
  +       !               integral approximations over the subintervals,
  +       !              work(limit*3+1), ..., work(limit*3)
  +       !               contain the error estimates.
  +       !
  +       implicit none
  +
  +       real ( kind = 8 ) abserr,bound,epsabs,epsrel,f,result,work
  +       integer ( kind = 4 ) ier,inf,iwork,last,lenw,limit,lvl,l1,l2,l3,neval
  +
  +       dimension iwork(limit),work(lenw)
  +
  +       external f
  +       !
  +       !  check validity of limit and lenw.
  +       !
  +       ier = 6
  +       neval = 0
  +       last = 0
  +       result = 0.0_real64
  +       abserr = 0.0_real64
  +       if(limit.lt.1.or.lenw.lt.limit*4) go to 10
  +       !
  +       !  prepare call for dqagie.
  +       !
  +       l1 = limit+1
  +       l2 = limit+l1
  +       l3 = limit+l2
  +
  +       call dqagie(f,bound,inf,epsabs,epsrel,limit,result,abserr, &
  +            neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
  +       !
  +       !  call error handler if necessary.
  +       !
  +       lvl = 0
  +10     if(ier.eq.6) lvl = 1
  +
  +       if(ier.ne.0) then
  +          call xerror('abnormal return from dqagi',26,ier,lvl)
  +       end if
  +
  +       return
  +     end subroutine dqagi
  +
  +     subroutine dqage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, &
  +          neval, ier, alist, blist, rlist, elist, iord, last )
  +
  +       !*****************************************************************************80
  +       !
  +       !! DQAGE estimates a definite integral.
  +       !
  +       !  Modified:
  +       !
  +       !    11 September 2015
  +       !
  +       !  Author:
  +       !
  +       !    Robert Piessens, Elise de Doncker
  +       !
  +       !***purpose  the routine calculates an approximation result to a given
  +       !      definite integral   i = integral of f over (a,b),
  +       !      hopefully satisfying following claim for accuracy
  +       !      abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
  +       !
  +       !  Parameters:
  +       !
  +       !   on entry
  +       !      f      - real(real64)
  +       !               function subprogram defining the integrand
  +       !               function f(x). the actual name for f needs to be
  +       !               declared e x t e r n a l in the driver program.
  +       !
  +       !      a      - real(real64)
  +       !               lower limit of integration
  +       !
  +       !      b      - real(real64)
  +       !               upper limit of integration
  +       !
  +       !      epsabs - real(real64)
  +       !               absolute accuracy requested
  +       !      epsrel - real(real64)
  +       !               relative accuracy requested
  +       !               if  epsabs.le.0
  +       !               and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
  +       !               the routine will end with ier = 6.
  +       !
  +       !      key    - integer(int32)
  +       !               key for choice of local integration rule
  +       !               a gauss-kronrod pair is used with
  +       !                    7 - 15 points if key.lt.2,
  +       !                   10 - 21 points if key = 2,
  +       !                   15 - 31 points if key = 3,
  +       !                   20 - 41 points if key = 4,
  +       !                   25 - 51 points if key = 5,
  +       !                   30 - 61 points if key.gt.5.
  +       !
  +       !      limit  - integer(int32)
  +       !               gives an upperbound on the number of subintervals
  +       !               in the partition of (a,b), limit.ge.1.
  +       !
  +       !   on return
  +       !      result - real(real64)
  +       !               approximation to the integral
  +       !
  +       !      abserr - real(real64)
  +       !               estimate of the modulus of the absolute error,
  +       !               which should equal or exceed abs(i-result)
  +       !
  +       !      neval  - integer(int32)
  +       !               number of integrand evaluations
  +       !
  +       !      ier    - integer(int32)
  +       !               ier = 0 normal and reliable termination of the
  +       !                       routine. it is assumed that the requested
  +       !                       accuracy has been achieved.
  +       !               ier.gt.0 abnormal termination of the routine
  +       !                       the estimates for result and error are
  +       !                       less reliable. it is assumed that the
  +       !                       requested accuracy has not been achieved.
  +       !      error messages
  +       !               ier = 1 maximum number of subdivisions allowed
  +       !                       has been achieved. one can allow more
  +       !                       subdivisions by increasing the value
  +       !                       of limit.
  +       !                       however, if this yields no improvement it
  +       !                       is rather advised to analyze the integrand
  +       !                       in order to determine the integration
  +       !                       difficulties. if the position of a local
  +       !                       difficulty can be determined(e.g.
  +       !                       singularity, discontinuity within the
  +       !                       interval) one will probably gain from
  +       !                       splitting up the interval at this point
  +       !                       and calling the integrator on the
  +       !                       subranges. if possible, an appropriate
  +       !                       special-purpose integrator should be used
  +       !                       which is designed for handling the type of
  +       !                       difficulty involved.
  +       !                   = 2 the occurrence of roundoff error is
  +       !                       detected, which prevents the requested
  +       !                       tolerance from being achieved.
  +       !                   = 3 extremely bad integrand behaviour occurs
  +       !                       at some points of the integration
  +       !                       interval.
  +       !                   = 6 the input is invalid, because
  +       !                       (epsabs.le.0 and
  +       !                        epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
  +       !                       result, abserr, neval, last, rlist(1) ,
  +       !                       elist(1) and iord(1) are set to zero.
  +       !                       alist(1) and blist(1) are set to a and b
  +       !                       respectively.
  +       !
  +       !      alist   - real(real64)
  +       !                vector of dimension at least limit, the first
  +       !                 last  elements of which are the left
  +       !                end points of the subintervals in the partition
  +       !                of the given integration range (a,b)
  +       !
  +       !      blist   - real(real64)
  +       !                vector of dimension at least limit, the first
  +       !                 last  elements of which are the right
  +       !                end points of the subintervals in the partition
  +       !                of the given integration range (a,b)
  +       !
  +       !      rlist   - real(real64)
  +       !                vector of dimension at least limit, the first
  +       !                 last  elements of which are the
  +       !                integral approximations on the subintervals
  +       !
  +       !      elist   - real(real64)
  +       !                vector of dimension at least limit, the first
  +       !                 last  elements of which are the moduli of the
  +       !                absolute error estimates on the subintervals
  +       !
  +       !      iord    - integer(int32)
  +       !                vector of dimension at least limit, the first k
  +       !                elements of which are pointers to the
  +       !                error estimates over the subintervals,
  +       !                such that elist(iord(1)), ...,
  +       !                elist(iord(k)) form a decreasing sequence,
  +       !                with k = last if last.le.(limit/2+2), and
  +       !                k = limit+1-last otherwise
  +       !
  +       !      last    - integer(int32)
  +       !                number of subintervals actually produced in the
  +       !                subdivision process
  +       !
  +       !  Local Parameters:
  +       !
  +       !     alist     - list of left end points of all subintervals
  +       !                 considered up to now
  +       !     blist     - list of right end points of all subintervals
  +       !                 considered up to now
  +       !     rlist(i)  - approximation to the integral over
  +       !                (alist(i),blist(i))
  +       !     elist(i)  - error estimate applying to rlist(i)
  +       !     maxerr    - pointer to the interval with largest
  +       !                 error estimate
  +       !     errmax    - elist(maxerr)
  +       !     area      - sum of the integrals over the subintervals
  +       !     errsum    - sum of the errors over the subintervals
  +       !     errbnd    - requested accuracy max(epsabs,epsrel*
  +       !                 abs(result))
  +       !     *****1    - variable for the left subinterval
  +       !     *****2    - variable for the right subinterval
  +       !     last      - index for subdivision
  +       !
  +       !
  +       !     machine dependent constants
  +       !
  +       !     epmach  is the largest relative spacing.
  +       !     uflow  is the smallest positive magnitude.
  +       !
  +       implicit none
  +
  +       real(real64) a,abserr,alist,area,area1,area12,area2,a1,a2,b, &
  +            blist,b1,b2,defabs,defab1,defab2,elist,epmach, &
  +            epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum, &
  +            resabs,result,rlist,uflow
  +       integer(int32) ier,iord,iroff1,iroff2,k,key,keyf,last,limit, &
  +            maxerr, nrmax, neval
  +       real(real64), external :: f
  +
  +       !    interface
  +       !       function f(x)
  +       !         use iso_fortran_env
  +       !         real(real64), intent(in) :: x
  +       !         real(real64) :: f
  +       !       end function f
  +       !    end interface
  +
  +       dimension alist(limit),blist(limit),elist(limit),iord(limit), &
  +            rlist(limit)
  +
  +       epmach = epsilon ( epmach )
  +       uflow = tiny ( uflow )
  +       !
  +       !  test on validity of parameters
  +       !
  +       ier = 0
  +       neval = 0
  +       last = 0
  +       result = 0.0_real64
  +       abserr = 0.0_real64
  +       alist(1) = a
  +       blist(1) = b
  +       rlist(1) = 0.0_real64
  +       elist(1) = 0.0_real64
  +       iord(1) = 0
  +
  +       if(epsabs.le.0.0_real64.and. &
  +            epsrel.lt. max ( 0.5e2_real64*epmach,0.5e-28_real64)) then
  +          ier = 6
  +          return
  +       end if
  +       !
  +       !  first approximation to the integral
  +       !
  +       keyf = key
  +       if(key.le.0) keyf = 1
  +       if(key.ge.7) keyf = 6
  +       neval = 0
  +       if(keyf.eq.1) call dqk15(f,a,b,result,abserr,defabs,resabs)
  +       if(keyf.eq.2) call dqk21(f,a,b,result,abserr,defabs,resabs)
  +       if(keyf.eq.3) call dqk31(f,a,b,result,abserr,defabs,resabs)
  +       if(keyf.eq.4) call dqk41(f,a,b,result,abserr,defabs,resabs)
  +       if(keyf.eq.5) call dqk51(f,a,b,result,abserr,defabs,resabs)
  +       if(keyf.eq.6) call dqk61(f,a,b,result,abserr,defabs,resabs)
  +       last = 1
  +       rlist(1) = result
  +       elist(1) = abserr
  +       iord(1) = 1
  +       !
  +       !  test on accuracy.
  +       !
  +       errbnd =  max ( epsabs, epsrel* abs ( result ) )
  +
  +       if(abserr.le.0.5e2_real64 * epmach * defabs .and. &
  +            abserr.gt.errbnd) then
  +          ier = 2
  +       end if
  +
  +       if(limit.eq.1) then
  +          ier = 1
  +       end if
  +
  +       if ( ier .ne. 0 .or. &
  +            (abserr .le. errbnd .and. abserr .ne. resabs ) .or. &
  +            abserr .eq. 0.0_real64 ) then
  +
  +          if(keyf.ne.1) then
  +             neval = (10*keyf+1)*(2*neval+1)
  +          else
  +             neval = 30*neval+15
  +          end if
  +
  +          return
  +
  +       end if
  +       !
  +       !  initialization
  +       !
  +       errmax = abserr
  +       maxerr = 1
  +       area = result
  +       errsum = abserr
  +       nrmax = 1
  +       iroff1 = 0
  +       iroff2 = 0
  +       !
  +       !  main do-loop
  +       !
  +       do last = 2, limit
  +          !
  +          !  bisect the subinterval with the largest error estimate.
  +          !
  +          a1 = alist(maxerr)
  +          b1 = 0.5_real64*(alist(maxerr)+blist(maxerr))
  +          a2 = b1
  +          b2 = blist(maxerr)
  +
  +          if(keyf.eq.1) call dqk15(f,a1,b1,area1,error1,resabs,defab1)
  +          if(keyf.eq.2) call dqk21(f,a1,b1,area1,error1,resabs,defab1)
  +          if(keyf.eq.3) call dqk31(f,a1,b1,area1,error1,resabs,defab1)
  +          if(keyf.eq.4) call dqk41(f,a1,b1,area1,error1,resabs,defab1)
  +          if(keyf.eq.5) call dqk51(f,a1,b1,area1,error1,resabs,defab1)
  +          if(keyf.eq.6) call dqk61(f,a1,b1,area1,error1,resabs,defab1)
  +
  +          if(keyf.eq.1) call dqk15(f,a2,b2,area2,error2,resabs,defab2)
  +          if(keyf.eq.2) call dqk21(f,a2,b2,area2,error2,resabs,defab2)
  +          if(keyf.eq.3) call dqk31(f,a2,b2,area2,error2,resabs,defab2)
  +          if(keyf.eq.4) call dqk41(f,a2,b2,area2,error2,resabs,defab2)
  +          if(keyf.eq.5) call dqk51(f,a2,b2,area2,error2,resabs,defab2)
  +          if(keyf.eq.6) call dqk61(f,a2,b2,area2,error2,resabs,defab2)
  +          !
  +          !  improve previous approximations to integral
  +          !  and error and test for accuracy.
  +          !
  +          neval = neval+1
  +          area12 = area1+area2
  +          erro12 = error1+error2
  +          errsum = errsum+erro12-errmax
  +          area = area+area12-rlist(maxerr)
  +
  +          if ( defab1 .ne. error1 .and. defab2 .ne. error2 ) then
  +
  +             if( abs ( rlist(maxerr)-area12).le.0.1e-4_real64* abs ( area12) &
  +                  .and. erro12.ge.0.99_real64*errmax) then
  +                iroff1 = iroff1+1
  +             end if
  +
  +             if(last.gt.10.and.erro12.gt.errmax) then
  +                iroff2 = iroff2+1
  +             end if
  +
  +          end if
  +
  +          rlist(maxerr) = area1
  +          rlist(last) = area2
  +          errbnd =  max ( epsabs,epsrel* abs ( area))
  +
  +          if ( errbnd .lt. errsum ) then
  +             !
  +             !  test for roundoff error and eventually set error flag.
  +             !
  +             if(iroff1.ge.6.or.iroff2.ge.20) then
  +                ier = 2
  +             end if
  +             !
  +             !  set error flag in the case that the number of subintervals
  +             !  equals limit.
  +             !
  +             if(last.eq.limit) then
  +                ier = 1
  +             end if
  +             !
  +             !  set error flag in the case of bad integrand behaviour
  +             !  at a point of the integration range.
  +             !
  +             if( max (  abs ( a1), abs ( b2)).le.(0.1e1_real64+0.1e3_real64* &
  +                  epmach)*( abs ( a2)+0.1e4_real64*uflow)) then
  +                ier = 3
  +             end if
  +
  +          end if
  +          !
  +          !  append the newly-created intervals to the list.
  +          !
  +          if(error2.le.error1) then
  +             alist(last) = a2
  +             blist(maxerr) = b1
  +             blist(last) = b2
  +             elist(maxerr) = error1
  +             elist(last) = error2
  +          else
  +             alist(maxerr) = a2
  +             alist(last) = a1
  +             blist(last) = b1
  +             rlist(maxerr) = area2
  +             rlist(last) = area1
  +             elist(maxerr) = error2
  +             elist(last) = error1
  +          end if
  +          !
  +          !  call dqpsrt to maintain the descending ordering
  +          !  in the list of error estimates and select the subinterval
  +          !  with the largest error estimate (to be bisected next).
  +          !
  +          call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
  +
  +          if(ier.ne.0.or.errsum.le.errbnd) then
  +             exit
  +          end if
  +
  +       end do
  +       !
  +       !  compute final result.
  +       !
  +       result = 0.0_real64
  +       do k=1,last
  +          result = result+rlist(k)
  +       end do
  +       abserr = errsum
   
          if(keyf.ne.1) then
             neval = (10*keyf+1)*(2*neval+1)
  
positive.code:     !    Robert Piessens, Elise de Doncker
      !
      !***purpose  the routine calculates an approximation result to a given
      !      definite integral   i = integral of f over (a,b),
      !      hopefully satisfying following claim for accuracy
      !      abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
      !
      !  Parameters:
      !
      !   on entry
      !      f      - real(real64)
      !               function subprogram defining the integrand
      !               function f(x). the actual name for f needs to be
      !               declared e x t e r n a l in the driver program.
      !
      !      a      - real(real64)
      !               lower limit of integration
      !
      !      b      - real(real64)
      !               upper limit of integration
      !
      !      epsabs - real(real64)
      !               absolute accuracy requested
      !      epsrel - real(real64)
      !               relative accuracy requested
      !               if  epsabs.le.0
      !               and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
      !               the routine will end with ier = 6.
      !
      !      key    - integer(int32)
      !               key for choice of local integration rule
      !               a gauss-kronrod pair is used with
      !                    7 - 15 points if key.lt.2,
      !                   10 - 21 points if key = 2,
      !                   15 - 31 points if key = 3,
      !                   20 - 41 points if key = 4,
      !                   25 - 51 points if key = 5,
      !                   30 - 61 points if key.gt.5.
      !
      !      limit  - integer(int32)
      !               gives an upperbound on the number of subintervals
      !               in the partition of (a,b), limit.ge.1.
      !
      !   on return
      !      result - real(real64)
      !               approximation to the integral
      !
      !      abserr - real(real64)
      !               estimate of the modulus of the absolute error,
      !               which should equal or exceed abs(i-result)
      !
      !      neval  - integer(int32)
      !               number of integrand evaluations
      !
      !      ier    - integer(int32)
      !               ier = 0 normal and reliable termination of the
      !                       routine. it is assumed that the requested
      !                       accuracy has been achieved.
      !               ier.gt.0 abnormal termination of the routine
      !                       the estimates for result and error are
      !                       less reliable. it is assumed that the
      !                       requested accuracy has not been achieved.
      !      error messages
      !               ier = 1 maximum number of subdivisions allowed
      !                       has been achieved. one can allow more
      !                       subdivisions by increasing the value
      !                       of limit.
      !                       however, if this yields no improvement it
      !                       is rather advised to analyze the integrand
      !                       in order to determine the integration
      !                       difficulties. if the position of a local
      !                       difficulty can be determined(e.g.
      !                       singularity, discontinuity within the
      !                       interval) one will probably gain from
      !                       splitting up the interval at this point
      !                       and calling the integrator on the
      !                       subranges. if possible, an appropriate
      !                       special-purpose integrator should be used
      !                       which is designed for handling the type of
      !                       difficulty involved.
      !                   = 2 the occurrence of roundoff error is
      !                       detected, which prevents the requested
      !                       tolerance from being achieved.
      !                   = 3 extremely bad integrand behaviour occurs
      !                       at some points of the integration
      !                       interval.
      !                   = 6 the input is invalid, because
      !                       (epsabs.le.0 and
      !                        epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
      !                       result, abserr, neval, last, rlist(1) ,
      !                       elist(1) and iord(1) are set to zero.
      !                       alist(1) and blist(1) are set to a and b
      !                       respectively.
      !
      !      alist   - real(real64)
      !                vector of dimension at least limit, the first
      !                 last  elements of which are the left
      !                end points of the subintervals in the partition
      !                of the given integration range (a,b)
      !
      !      blist   - real(real64)
      !                vector of dimension at least limit, the first
      !                 last  elements of which are the right
      !                end points of the subintervals in the partition
      !                of the given integration range (a,b)
      !
      !      rlist   - real(real64)
      !                vector of dimension at least limit, the first
      !                 last  elements of which are the
      !                integral approximations on the subintervals
      !
      !      elist   - real(real64)
      !                vector of dimension at least limit, the first
      !                 last  elements of which are the moduli of the
      !                absolute error estimates on the subintervals
      !
      !      iord    - integer(int32)
      !                vector of dimension at least limit, the first k
      !                elements of which are pointers to the
      !                error estimates over the subintervals,
      !                such that elist(iord(1)), ...,
      !                elist(iord(k)) form a decreasing sequence,
      !                with k = last if last.le.(limit/2+2), and
      !                k = limit+1-last otherwise
      !
      !      last    - integer(int32)
      !                number of subintervals actually produced in the
      !                subdivision process
      !
      !  Local Parameters:
      !
      !     alist     - list of left end points of all subintervals
      !                 considered up to now
      !     blist     - list of right end points of all subintervals
      !                 considered up to now
      !     rlist(i)  - approximation to the integral over
      !                (alist(i),blist(i))
      !     elist(i)  - error estimate applying to rlist(i)
      !     maxerr    - pointer to the interval with largest
      !                 error estimate
      !     errmax    - elist(maxerr)
      !     area      - sum of the integrals over the subintervals
      !     errsum    - sum of the errors over the subintervals
      !     errbnd    - requested accuracy max(epsabs,epsrel*
  
negative.code:   implicit none
    ! 15 June 2016
    ! quadrature module
    ! routines from QUADPACK, fortran90 version:
    ! https://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.f90
    ! kind types parameters from the iso_fortran_env intrinsic module
    ! Author: Simone Marsili
  
  contains
  
    subroutine dqagie ( f, bound, inf, epsabs, epsrel, limit, result, abserr, &
         neval, ier, alist, blist, rlist, elist, iord, last )
  
      !*****************************************************************************80
      !
      !! DQAGIE estimates an integral over a semi-infinite or infinite interval.
      !
      !  Modified:
      !
  