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:     !     *****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
      !
  
negative.code:     !    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:
      !
      !      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    - real ( kind = 8 )
      !               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.
      !
      !      limit  - integer ( kind = 4 )
      !               gives an upper bound on the number of subintervals
      !               in the partition of (a,b), limit.ge.1
      !
      !   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),
      !                       result, abserr, neval, last, rlist(1),
      !                       elist(1) and iord(1) are set to zero.
      !                       alist(1) and blist(1) are set to 0
      !                       and 1 respectively.
      !
      !      alist  - real ( kind = 8 )
      !               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 transformed integration range (0,1).
      !
      !      blist  - real ( kind = 8 )
      !               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 transformed integration range (0,1).
      !
      !      rlist  - real ( kind = 8 )
      !               vector of dimension at least limit, the first
      !                last  elements of which are the integral
      !               approximations on the subintervals
      !
      !      elist  - real ( kind = 8 )
      !               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 ( kind = 4 )
      !               vector of dimension 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 ( kind = 4 )
      !               number of subintervals actually produced
      !               in the subdivision process
      !
      !  Local Parameters:
      !
      !      the dimension of rlist2 is determined by the value of
      !      limexp in routine dqelg.
      !
      !     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))
      !     rlist2    - array of dimension at least (limexp+2),
      !                 containing the part of the epsilon table
      !                 wich is still needed for further computations
      !     elist(i)  - error estimate applying to rlist(i)
      !     maxerr    - pointer to the interval with largest error
      !                 estimate
      !     errmax    - elist(maxerr)
      !     erlast    - error on the interval currently subdivided
      !                 (before that subdivision has taken place)
      !     area      - sum of the integrals over the subintervals
      !     errsum    - sum of the errors over the subintervals
      !     errbnd    - requested accuracy max(epsabs,epsrel*
  