anchor.diff:    implicit none
     ! 15 June 2016
     ! quadrature module
  -  ! routines from QUADPACK, fortran90 version: 
  +  ! 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 dqage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, &
  +  subroutine dqagie ( f, bound, inf, epsabs, epsrel, limit, result, abserr, &
          neval, ier, alist, blist, rlist, elist, iord, last )
   
       !*****************************************************************************80
       !
  -    !! DQAGE estimates a definite integral.
  +    !! DQAGIE estimates an integral over a semi-infinite or infinite interval.
       !
       !  Modified:
       !
  
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:     !     *****1    - variable for the left subinterval
      !     *****2    - variable for the right subinterval
      !     last      - index for subdivision
      !     nres      - number of calls to the extrapolation routine
      !     numrl2    - number of elements currently in rlist2. if an
      !                 appropriate approximation to the compounded
      !                 integral has been obtained, it is put in
      !                 rlist2(numrl2) after numrl2 has been increased
      !                 by one.
      !     small     - length of the smallest interval considered up
      !                 to now, multiplied by 1.5
      !     erlarg    - sum of the errors over the intervals larger
      !                 than the smallest interval considered up to now
      !     extrap    - logical variable denoting that the routine
      !                 is attempting to perform extrapolation. i.e.
      !                 before subdividing the smallest interval we
      !                 try to decrease the value of erlarg.
      !     noext     - logical variable denoting that extrapolation
      !                 is no longer allowed (true-value)
      !
      !      machine dependent constants
      !
      !     epmach is the largest relative spacing.
      !     uflow is the smallest positive magnitude.
      !     oflow is the largest positive magnitude.
      !
      implicit none
  
      real (real64) abseps,abserr,alist,area,area1,area12,area2,a1, &
           a2,blist,boun,bound,b1,b2,correc,defabs,defab1,defab2, &
           dres,elist,epmach,epsabs,epsrel,erlarg,erlast, &
           errbnd,errmax,error1,error2,erro12,errsum,ertest,f,oflow,resabs, &
           reseps,result,res3la,rlist,rlist2,small,uflow
      integer (int32) id,ier,ierro,inf,iord,iroff1,iroff2, &
           iroff3,jupbnd,k,ksgn, &
           ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2
      logical extrap,noext
      dimension alist(limit),blist(limit),elist(limit),iord(limit), &
           res3la(3),rlist(limit),rlist2(52)
  
      external f
  
      epmach = epsilon ( epmach )
      !
      !  test on validity of parameters
      !
  