Line data    Source code

       1             : !! ******************************************************************
       2             : !! The routines listed in this file "adgaquad_mod.F90" are performing
       3             : !! Numerical Integrations using some kind of 
       4             : !! adaptive Gauss quadrature.
       5             : !! They are taken from the Internet (http://www.netlib.org)
       6             : !! and parts of different software packages / libraries.
       7             : !! ******************************************************************
       8             : !! For any restrictions on the use of the routines, please see
       9             : !! the original web site. 
      10             : !! ******************************************************************
      11             : !! Changes: calls to error handler 'xerror()' replaced by
      12             : !!          WRITE(7,*) - statements.
      13             : !! ******************************************************************
      14             : !! list of routines and the libraries they are taken from:
      15             : !!  dqag                calling routine, bounded integration interval
      16             : !!                      QUADPACK; calls: dqage
      17             : !!  dqage               the integration routine, bounded interval
      18             : !!                      QUADPACK; calls: sd1mach,dqk15,dqk21,dqk31,
      19             : !!                                      dqk41,dqk51,dqk61,dqpsrt
      20             : !!  dqagi               calling routine, unbounded (semi-infinite or 
      21             : !!                      infinite) integration interval
      22             : !!                      QUADPACK; calls: dqagie
      23             : !!  dqagie              the integration routine, unbounded interval
      24             : !!                      QUADPACK; calls: sd1mach,dqelg,dqk15i,dqpsrt
      25             : !! ------------------------------------------------------------------
      26             : !!  dqk15               QUADPACK; calls: sd1mach
      27             : !!  dqk21               QUADPACK; calls: sd1mach
      28             : !!  dqk31               QUADPACK; calls: sd1mach
      29             : !!  dqk41               QUADPACK; calls: sd1mach
      30             : !!  dqk51               QUADPACK; calls: sd1mach
      31             : !!  dqk61               QUADPACK; calls: sd1mach
      32             : !!  dqpsrt              QUADPACK; calls: none
      33             : !!  dqk15i              QUADPACK; calls: sd1mach
      34             : !!  dqelg               QUADPACK; calls: sd1mach
      35             : !! ------------------------------------------------------------------
      36             : !!  xerror              Error handling routine
      37             : !!                      ALLIANT (/quad); calls: xerrwv
      38             : !!  xerrwv              Error handling routine
      39             : !!                      SODEPACK; calls: none
      40             : !!  d1mach              determine machine parameters (accuracies)
      41             : !!                      BLAS; calls: none
      42             : !! ------------------------------------------------------------------
      43             : 
      44             : module adgaquad_mod
      45             : 
      46             :   use carma_precision_mod
      47             :   use adgaquad_types_mod
      48             : 
      49             :   implicit none
      50             :   
      51             :   private
      52             :  
      53             :   public :: dqag
      54             :   public :: dqage
      55             :   public :: dqagi
      56             :   public :: dqagie
      57             :   
      58             :   contains 
      59             : 
      60             :   !!***begin prologue  dqag
      61             :   !!***date written   800101   (yymmdd)
      62             :   !!***revision date  130319   (yymmdd)  
      63             :   !!***category no.  h2a1a1
      64             :   !!***keywords  automatic integrator, general-purpose,
      65             :   !!             integrand examinator, globally adaptive,
      66             :   !!             gauss-kronrod
      67             :   !!***author  piessens,robert,appl. math. & progr. div - k.u.leuven
      68             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
      69             :   !!***purpose  the routine calculates an approximation result to a given
      70             :   !!            definite integral i = integral of f over (a,b),
      71             :   !!            hopefully satisfying following claim for accuracy
      72             :   !!            abs(i-result)le.max(epsabs,epsrel*abs(i)).
      73             :   !!***description
      74             :   !!
      75             :   !!        computation of a definite integral
      76             :   !!        standard fortran subroutine
      77             :   !!        double precision version
      78             :   !!
      79             :   !!            fx     - double precision
      80             :   !!                     function subprogam defining the integrand
      81             :   !!                     function f(x). the actual name for f needs to be
      82             :   !!                     declared e x t e r n a l in the driver program.
      83             :   !!
      84             :   !!            fx_vars- structure containing variables need for integration
      85             :   !!                     specific to fractal meanfield scattering code
      86             :   !!
      87             :   !!            a      - double precision
      88             :   !!                     lower limit of integration
      89             :   !!
      90             :   !!            b      - double precision
      91             :   !!                     upper limit of integration
      92             :   !!
      93             :   !!            epsabs - double precision
      94             :   !!                     absolute accoracy requested
      95             :   !!            epsrel - double precision
      96             :   !!                     relative accuracy requested
      97             :   !!                     if  epsabs.le.0
      98             :   !!                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
      99             :   !!                     the routine will end with ier = 6.
     100             :   !!
     101             :   !!            key    - integer
     102             :   !!                     key for choice of local integration rule
     103             :   !!                     a gauss-kronrod pair is used with
     104             :   !!                       7 - 15 points if key.lt.2,
     105             :   !!                      10 - 21 points if key = 2,
     106             :   !!                      15 - 31 points if key = 3,
     107             :   !!                      20 - 41 points if key = 4,
     108             :   !!                      25 - 51 points if key = 5,
     109             :   !!                      30 - 61 points if key.gt.5.
     110             :   !!
     111             :   !!         on return
     112             :   !!            result - double precision
     113             :   !!                     approximation to the integral
     114             :   !!
     115             :   !!            abserr - double precision
     116             :   !!                     estimate of the modulus of the absolute error,
     117             :   !!                     which should equal or exceed abs(i-result)
     118             :   !!
     119             :   !!            neval  - integer
     120             :   !!                     number of integrand evaluations
     121             :   !!
     122             :   !!            ier    - integer
     123             :   !!                     ier = 0 normal and reliable termination of the
     124             :   !!                             routine. it is assumed that the requested
     125             :   !!                             accuracy has been achieved.
     126             :   !!                     ier.gt.0 abnormal termination of the routine
     127             :   !!                             the estimates for result and error are
     128             :   !!                             less reliable. it is assumed that the
     129             :   !!                             requested accuracy has not been achieved.
     130             :   !!                      error messages
     131             :   !!                     ier = 1 maximum number of subdivisions allowed
     132             :   !!                             has been achieved. one can allow more
     133             :   !!                             subdivisions by increasing the value of
     134             :   !!                             limit (and taking the according dimension
     135             :   !!                             adjustments into account). however, if
     136             :   !!                             this yield no improvement it is advised
     137             :   !!                             to analyze the integrand in order to
     138             :   !!                             determine the integration difficulaties.
     139             :   !!                             if the position of a local difficulty can
     140             :   !!                             be determined (i.e.singularity,
     141             :   !!                             discontinuity within the interval) one
     142             :   !!                             will probably gain from splitting up the
     143             :   !!                             interval at this point and calling the
     144             :   !!                             integrator on the subranges. if possible,
     145             :   !!                             an appropriate special-purpose integrator
     146             :   !!                             should be used which is designed for
     147             :   !!                             handling the type of difficulty involved.
     148             :   !!                         = 2 the occurrence of roundoff error is
     149             :   !!                             detected, which prevents the requested
     150             :   !!                             tolerance from being achieved.
     151             :   !!                         = 3 extremely bad integrand behaviour occurs
     152             :   !!                             at some points of the integration
     153             :   !!                             interval.
     154             :   !!                         = 6 the input is invalid, because
     155             :   !!                             (epsabs.le.0 and
     156             :   !!                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
     157             :   !!                             or limit.lt.1 or lenw.lt.limit*4.
     158             :   !!                             result, abserr, neval, last are set
     159             :   !!                             to zero.
     160             :   !!                             except when lenw is invalid, iwork(1),
     161             :   !!                             work(limit*2+1) and work(limit*3+1) are
     162             :   !!                             set to zero, work(1) is set to a and
     163             :   !!                             work(limit+1) to b.
     164             :   !!                         = 9 failure in sd1mach determining machine parameters
     165             :   !!
     166             :   !!         dimensioning parameters
     167             :   !!            limit - integer
     168             :   !!                    dimensioning parameter for iwork
     169             :   !!                    limit determines the maximum number of subintervals
     170             :   !!                    in the partition of the given integration interval
     171             :   !!                    (a,b), limit.ge.1.
     172             :   !!                    if limit.lt.1, the routine will end with ier = 6.
     173             :   !!
     174             :   !!            lenw  - integer
     175             :   !!                    dimensioning parameter for work
     176             :   !!                    lenw must be at least limit*4.
     177             :   !!                    if lenw.lt.limit*4, the routine will end with
     178             :   !!                    ier = 6.
     179             :   !!
     180             :   !!            last  - integer
     181             :   !!                    on return, last equals the number of subintervals
     182             :   !!                    produced in the subdiviosion process, which
     183             :   !!                    determines the number of significant elements
     184             :   !!                    actually in the work arrays.
     185             :   !!
     186             :   !!         work arrays
     187             :   !!            iwork - integer
     188             :   !!                    vector of dimension at least limit, the first k
     189             :   !!                    elements of which contain pointers to the error
     190             :   !!                    estimates over the subintervals, such that
     191             :   !!                    work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
     192             :   !!                    form a decreasing sequence with k = last if
     193             :   !!                    last.le.(limit/2+2), and k = limit+1-last otherwise
     194             :   !!
     195             :   !!            work  - double precision
     196             :   !!                    vector of dimension at least lenw
     197             :   !!                    on return
     198             :   !!                    work(1), ..., work(last) contain the left end
     199             :   !!                    points of the subintervals in the partition of
     200             :   !!                     (a,b),
     201             :   !!                    work(limit+1), ..., work(limit+last) contain the
     202             :   !!                     right end points,
     203             :   !!                    work(limit*2+1), ..., work(limit*2+last) contain
     204             :   !!                     the integral approximations over the subintervals,
     205             :   !!                    work(limit*3+1), ..., work(limit*3+last) contain
     206             :   !!                     the error estimates.
     207             :   !!
     208             :   !!***references  (none)
     209             :   !!***routines called  dqage,xerror
     210             :   !!***end prologue  dqag
     211           0 :   subroutine dqag(fx,fx_vars,a,b,epsabs,epsrel,key,result,abserr,neval,ier, &
     212           0 :                   limit,lenw,last,iwork,work)
     213             : 
     214             :     ! Arguments
     215             :     interface
     216             :       function fx(centr, vars)
     217             :         use carma_precision_mod, only : f
     218             :         use adgaquad_types_mod
     219             :         real(kind=f), intent(in) :: centr
     220             :         type(adgaquad_vars_type), intent(inout) :: vars
     221             :         real(kind=f)             :: fx
     222             :       end function fx
     223             :     end interface
     224             :     type(adgaquad_vars_type) :: fx_vars
     225             :     real(kind=f) :: a
     226             :     real(kind=f) :: b
     227             :     real(kind=f) :: epsabs
     228             :     real(kind=f) :: epsrel
     229             :     integer :: key
     230             :     real(kind=f) :: result
     231             :     real(kind=f) :: abserr
     232             :     integer :: neval
     233             :     integer :: ier
     234             :     integer :: limit
     235             :     integer :: lenw
     236             :     integer :: last
     237             :     integer :: iwork(limit)
     238             :     real(kind=f) :: work(lenw)
     239             : 
     240             :     ! Local declarations
     241             :     integer :: lvl,l1,l2,l3
     242             : 
     243             :     !  check validity of lenw.
     244             :     !
     245             :     !***first executable statement  dqag
     246           0 :     ier = 6
     247           0 :     neval = 0
     248           0 :     last = 0
     249           0 :     result = 0.0_f
     250           0 :     abserr = 0.0_f
     251           0 :     if(limit.lt.1.or.lenw.lt.limit*4) go to 10
     252             : 
     253             :     !  prepare call for dqage.
     254             : 
     255           0 :     l1 = limit+1
     256           0 :     l2 = limit+l1
     257           0 :     l3 = limit+l2
     258             : 
     259             :     call dqage(fx,fx_vars,a,b,epsabs,epsrel,key,limit,result,abserr,neval, &
     260           0 :                ier,work(1),work(l1),work(l2),work(l3),iwork,last)
     261             :  
     262             :     !  call error handler if necessary.
     263             :  
     264           0 :     lvl = 0
     265           0 : 10  if(ier.eq.6) lvl = 1
     266           0 :     if(ier.ne.0) then
     267           0 :       write(*,*) "ERROR: abnormal return from dqag"
     268           0 :       write(*,*) "       ifail=",ier,"  level=",lvl
     269             :     endif
     270           0 :     return
     271             :   end subroutine dqag
     272             : 
     273             : 
     274             : 
     275             :   !!***begin prologue  dqage
     276             :   !!***date written   800101   (yymmdd)
     277             :   !!***revision date  130319   (yymmdd)
     278             :   !!***category no.  h2a1a1
     279             :   !!***keywords  automatic integrator, general-purpose,
     280             :   !!             integrand examinator, globally adaptive,
     281             :   !!             gauss-kronrod
     282             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
     283             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
     284             :   !!***purpose  the routine calculates an approximation result to a given
     285             :   !!            definite integral   i = integral of f over (a,b),
     286             :   !!            hopefully satisfying following claim for accuracy
     287             :   !!            abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
     288             :   !!***description
     289             :   !!
     290             :   !!        computation of a definite integral
     291             :   !!        standard fortran subroutine
     292             :   !!        double precision version
     293             :   !!
     294             :   !!        parameters
     295             :   !!         on entry
     296             :   !!            fx     - double precision
     297             :   !!                     function subprogram defining the integrand
     298             :   !!                     function f(x). the actual name for f needs to be
     299             :   !!                     declared e x t e r n a l in the driver program.
     300             :   !!
     301             :   !!            fx_vars- structure containing variables need for integration
     302             :   !!                     specific to fractal meanfield scattering code
     303             :   !!
     304             :   !!            a      - double precision
     305             :   !!                     lower limit of integration
     306             :   !!
     307             :   !!            b      - double precision
     308             :   !!                     upper limit of integration
     309             :   !!
     310             :   !!            epsabs - double precision
     311             :   !!                     absolute accuracy requested
     312             :   !!            epsrel - double precision
     313             :   !!                     relative accuracy requested
     314             :   !!                     if  epsabs.le.0
     315             :   !!                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
     316             :   !!                     the routine will end with ier = 6.
     317             :   !!
     318             :   !!            key    - integer
     319             :   !!                     key for choice of local integration rule
     320             :   !!                     a gauss-kronrod pair is used with
     321             :   !!                          7 - 15 points if key.lt.2,
     322             :   !!                         10 - 21 points if key = 2,
     323             :   !!                         15 - 31 points if key = 3,
     324             :   !!                         20 - 41 points if key = 4,
     325             :   !!                         25 - 51 points if key = 5,
     326             :   !!                         30 - 61 points if key.gt.5.
     327             :   !!
     328             :   !!            limit  - integer
     329             :   !!                     gives an upperbound on the number of subintervals
     330             :   !!                     in the partition of (a,b), limit.ge.1.
     331             :   !!
     332             :   !!         on return
     333             :   !!            result - double precision
     334             :   !!                     approximation to the integral
     335             :   !!
     336             :   !!            abserr - double precision
     337             :   !!                     estimate of the modulus of the absolute error,
     338             :   !!                     which should equal or exceed abs(i-result)
     339             :   !!
     340             :   !!            neval  - integer
     341             :   !!                     number of integrand evaluations
     342             :   !!
     343             :   !!            ier    - integer
     344             :   !!                     ier = 0 normal and reliable termination of the
     345             :   !!                             routine. it is assumed that the requested
     346             :   !!                             accuracy has been achieved.
     347             :   !!                     ier.gt.0 abnormal termination of the routine
     348             :   !!                             the estimates for result and error are
     349             :   !!                             less reliable. it is assumed that the
     350             :   !!                             requested accuracy has not been achieved.
     351             :   !!            error messages
     352             :   !!                     ier = 1 maximum number of subdivisions allowed
     353             :   !!                             has been achieved. one can allow more
     354             :   !!                             subdivisions by increasing the value
     355             :   !!                             of limit.
     356             :   !!                             however, if this yields no improvement it
     357             :   !!                             is rather advised to analyze the integrand
     358             :   !!                             in order to determine the integration
     359             :   !!                             difficulties. if the position of a local
     360             :   !!                             difficulty can be determined(e.g.
     361             :   !!                             singularity, discontinuity within the
     362             :   !!                             interval) one will probably gain from
     363             :   !!                             splitting up the interval at this point
     364             :   !!                             and calling the integrator on the
     365             :   !!                             subranges. if possible, an appropriate
     366             :   !!                             special-purpose integrator should be used
     367             :   !!                             which is designed for handling the type of
     368             :   !!                             difficulty involved.
     369             :   !!                         = 2 the occurrence of roundoff error is
     370             :   !!                             detected, which prevents the requested
     371             :   !!                             tolerance from being achieved.
     372             :   !!                         = 3 extremely bad integrand behaviour occurs
     373             :   !!                             at some points of the integration
     374             :   !!                             interval.
     375             :   !!                         = 6 the input is invalid, because
     376             :   !!                             (epsabs.le.0 and
     377             :   !!                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
     378             :   !!                             result, abserr, neval, last, rlist(1) ,
     379             :   !!                             elist(1) and iord(1) are set to zero.
     380             :   !!                             alist(1) and blist(1) are set to a and b
     381             :   !!                             respectively.
     382             :   !!                         = 9 failure in sd1mach determining machine parameters
     383             :   !!
     384             :   !!            alist   - double precision
     385             :   !!                      vector of dimension at least limit, the first
     386             :   !!                       last  elements of which are the left
     387             :   !!                      end points of the subintervals in the partition
     388             :   !!                      of the given integration range (a,b)
     389             :   !!
     390             :   !!            blist   - double precision
     391             :   !!                      vector of dimension at least limit, the first
     392             :   !!                       last  elements of which are the right
     393             :   !!                      end points of the subintervals in the partition
     394             :   !!                      of the given integration range (a,b)
     395             :   !!
     396             :   !!            rlist   - double precision
     397             :   !!                      vector of dimension at least limit, the first
     398             :   !!                       last  elements of which are the
     399             :   !!                      integral approximations on the subintervals
     400             :   !!
     401             :   !!            elist   - double precision
     402             :   !!                      vector of dimension at least limit, the first
     403             :   !!                       last  elements of which are the moduli of the
     404             :   !!                       absolute error estimates on the subintervals
     405             :   !!
     406             :   !!            iord    - integer
     407             :   !!                      vector of dimension at least limit, the first k
     408             :   !!                      elements of which are pointers to the
     409             :   !!                      error estimates over the subintervals,
     410             :   !!                      such that elist(iord(1)), ...,
     411             :   !!                      elist(iord(k)) form a decreasing sequence,
     412             :   !!                      with k = last if last.le.(limit/2+2), and
     413             :   !!                      k = limit+1-last otherwise
     414             :   !!
     415             :   !!            last    - integer
     416             :   !!                      number of subintervals actually produced in the
     417             :   !!                      subdivision process
     418             :   !!
     419             :   !!***references  (none)
     420             :   !!***routines called  sd1mach,dqk15,dqk21,dqk31,
     421             :   !!                    dqk41,dqk51,dqk61,dqpsrt
     422             :   !!***end prologue  dqage
     423           0 :   subroutine dqage(fx,fx_vars,a,b,epsabs,epsrel,key,limit,result,abserr, &
     424           0 :                    neval,ier,alist,blist,rlist,elist,iord,last)
     425             :   
     426             :     ! Arguments
     427             :     interface
     428             :       function fx(centr, vars)
     429             :         use carma_precision_mod, only : f
     430             :         use adgaquad_types_mod
     431             :         real(kind=f), intent(in) :: centr
     432             :         type(adgaquad_vars_type), intent(inout) :: vars
     433             :         real(kind=f)             :: fx
     434             :       end function fx
     435             :     end interface
     436             :     type(adgaquad_vars_type) :: fx_vars
     437             :     real(kind=f) :: a
     438             :     real(kind=f) :: b
     439             :     real(kind=f) :: epsabs
     440             :     real(kind=f) :: epsrel
     441             :     integer :: limit
     442             :     integer :: key
     443             :     real(kind=f) :: result
     444             :     real(kind=f) :: abserr
     445             :     integer :: neval
     446             :     integer :: ier
     447             :     real(kind=f) :: alist(limit)
     448             :     real(kind=f) :: blist(limit)
     449             :     real(kind=f) :: rlist(limit)
     450             :     real(kind=f) :: elist(limit)
     451             :     integer :: iord(limit)
     452             :     integer :: last
     453             : 
     454             :     ! Local declarations
     455             :     real(kind=f) :: area, area1, area12, area2, a1, a2
     456             :     real(kind=f) :: b1, b2, dabs, defabs, defab1, defab2, dmax1, epmach   
     457             :     real(kind=f) :: errbnd,errmax,error1,error2,erro12,errsum,resabs,uflow
     458             :     integer :: iroff1,iroff2,k,keyf,maxerr,nrmax
     459             : 
     460             : 
     461             :     !            list of major variables
     462             :     !            -----------------------
     463             :     !
     464             :     !           alist     - list of left end points of all subintervals
     465             :     !                       considered up to now
     466             :     !           blist     - list of right end points of all subintervals
     467             :     !                       considered up to now
     468             :     !           rlist(i)  - approximation to the integral over
     469             :     !                      (alist(i),blist(i))
     470             :     !           elist(i)  - error estimate applying to rlist(i)
     471             :     !           maxerr    - pointer to the interval with largest
     472             :     !                       error estimate
     473             :     !           errmax    - elist(maxerr)
     474             :     !           area      - sum of the integrals over the subintervals
     475             :     !           errsum    - sum of the errors over the subintervals
     476             :     !           errbnd    - requested accuracy max(epsabs,epsrel*
     477             :     !                       abs(result))
     478             :     !           *****1    - variable for the left subinterval
     479             :     !           *****2    - variable for the right subinterval
     480             :     !           last      - index for subdivision
     481             :     !
     482             :     !
     483             :     !           machine dependent constants
     484             :     !           ---------------------------
     485             :     !
     486             :     !           epmach  is the largest relative spacing.
     487             :     !           uflow  is the smallest positive magnitude.
     488             :     !
     489             :     !***first executable statement  dqage
     490           0 :     call sd1mach(4,epmach,ier)
     491           0 :     if(ier.eq.9) return
     492           0 :     call sd1mach(1,uflow,ier)
     493           0 :     if(ier.eq.9) return
     494             :  
     495             :     !epmach = d1mach(4)
     496             :     !uflow = d1mach(1)
     497             : 
     498             :     !           test on validity of parameters
     499             :     !           ------------------------------
     500             :     !
     501           0 :     ier = 0
     502           0 :     neval = 0
     503           0 :     last = 0
     504           0 :     result = 0.0_f
     505           0 :     abserr = 0.0_f
     506           0 :     alist(1) = a
     507           0 :     blist(1) = b
     508           0 :     rlist(1) = 0.0_f
     509           0 :     elist(1) = 0.0_f
     510           0 :     iord(1) = 0
     511           0 :     if(epsabs.le.0.0_f.and.epsrel.lt.dmax1(0.5e2_f*epmach,0.5e-28_f)) ier = 6
     512           0 :     if(ier.eq.6) go to 999
     513             : 
     514             :     !           first approximation to the integral
     515             :     !           -----------------------------------
     516             :     !
     517           0 :     keyf = key
     518           0 :     if(key.le.0) keyf = 1
     519           0 :     if(key.ge.7) keyf = 6
     520           0 :     neval = 0
     521           0 :     if(keyf.eq.1) call dqk15(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     522           0 :     if(ier.eq.9) return
     523           0 :     if(keyf.eq.2) call dqk21(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     524           0 :     if(ier.eq.9) return
     525           0 :     if(keyf.eq.3) call dqk31(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     526           0 :     if(ier.eq.9) return
     527           0 :     if(keyf.eq.4) call dqk41(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     528           0 :     if(ier.eq.9) return
     529           0 :     if(keyf.eq.5) call dqk51(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     530           0 :     if(ier.eq.9) return
     531           0 :     if(keyf.eq.6) call dqk61(fx,fx_vars,a,b,result,abserr,defabs,resabs,ier)
     532           0 :     if(ier.eq.9) return
     533           0 :     last = 1
     534           0 :     rlist(1) = result
     535           0 :     elist(1) = abserr
     536           0 :     iord(1) = 1
     537             :     !
     538             :     !           test on accuracy.
     539             :     !
     540           0 :     errbnd = dmax1(epsabs,epsrel*dabs(result))
     541           0 :     if(abserr.le.0.5e2_f*epmach*defabs.and.abserr.gt.errbnd) ier = 2
     542           0 :     if(limit.eq.1) ier = 1
     543           0 :     if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or.abserr.eq.0.0_f) go to 60
     544             :     !
     545             :     !           initialization
     546             :     !           --------------
     547             :     !
     548           0 :     errmax = abserr
     549           0 :     maxerr = 1
     550           0 :     area = result
     551           0 :     errsum = abserr
     552           0 :     nrmax = 1
     553           0 :     iroff1 = 0
     554           0 :     iroff2 = 0
     555             :     !
     556             :     !           main do-loop
     557             :     !           ------------
     558             :     !
     559           0 :     do last = 2,limit
     560             :       !
     561             :       !           bisect the subinterval with the largest error estimate.
     562             :       !
     563           0 :       a1 = alist(maxerr)
     564           0 :       b1 = 0.5_f*(alist(maxerr)+blist(maxerr))
     565           0 :       a2 = b1
     566           0 :       b2 = blist(maxerr)
     567           0 :       if(keyf.eq.1) call dqk15(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     568           0 :       if(ier.eq.9) return
     569           0 :       if(keyf.eq.2) call dqk21(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     570           0 :       if(ier.eq.9) return
     571           0 :       if(keyf.eq.3) call dqk31(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     572           0 :       if(ier.eq.9) return
     573           0 :       if(keyf.eq.4) call dqk41(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     574           0 :       if(ier.eq.9) return
     575           0 :       if(keyf.eq.5) call dqk51(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     576           0 :       if(ier.eq.9) return
     577           0 :       if(keyf.eq.6) call dqk61(fx,fx_vars,a1,b1,area1,error1,resabs,defab1,ier)
     578           0 :       if(ier.eq.9) return
     579           0 :       if(keyf.eq.1) call dqk15(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     580           0 :       if(ier.eq.9) return
     581           0 :       if(keyf.eq.2) call dqk21(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     582           0 :       if(ier.eq.9) return
     583           0 :       if(keyf.eq.3) call dqk31(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     584           0 :       if(ier.eq.9) return
     585           0 :       if(keyf.eq.4) call dqk41(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     586           0 :       if(ier.eq.9) return
     587           0 :       if(keyf.eq.5) call dqk51(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     588           0 :       if(ier.eq.9) return
     589           0 :       if(keyf.eq.6) call dqk61(fx,fx_vars,a2,b2,area2,error2,resabs,defab2,ier)
     590           0 :       if(ier.eq.9) return
     591             :       !           improve previous approximations to integral
     592             :       !           and error and test for accuracy.
     593             :       !
     594             :       !       neval = neval+1
     595           0 :       area12 = area1+area2
     596           0 :       erro12 = error1+error2
     597           0 :       errsum = errsum+erro12-errmax
     598           0 :       area = area+area12-rlist(maxerr)
     599           0 :       if(defab1.eq.error1.or.defab2.eq.error2) go to 5
     600           0 :       if(dabs(rlist(maxerr)-area12).le.0.1e-4_f*dabs(area12).and.erro12.ge.0.99_f*errmax) iroff1 = iroff1+1
     601           0 :       if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1
     602           0 :     5 rlist(maxerr) = area1
     603           0 :       rlist(last) = area2
     604           0 :       errbnd = dmax1(epsabs,epsrel*dabs(area))
     605           0 :       if(errsum.le.errbnd) go to 8
     606             :       !
     607             :       !           test for roundoff error and eventually set error flag.
     608             :       !
     609           0 :       if(iroff1.ge.6.or.iroff2.ge.20) ier = 2
     610             :       !
     611             :       !           set error flag in the case that the number of subintervals
     612             :       !           equals limit.
     613             :       !
     614           0 :       if(last.eq.limit) ier = 1
     615             :       !
     616             :       !           set error flag in the case of bad integrand behaviour
     617             :       !           at a point of the integration range.
     618             :       !
     619           0 :       if(dmax1(dabs(a1),dabs(b2)).le.(0.1e1_f+0.1e3_f*epmach)*(dabs(a2)+0.1e4_f*uflow)) ier = 3
     620             :       !
     621             :       !           append the newly-created intervals to the list.
     622             :       !
     623           0 :     8 if(error2.gt.error1) go to 10
     624           0 :       alist(last) = a2
     625           0 :       blist(maxerr) = b1
     626           0 :       blist(last) = b2
     627           0 :       elist(maxerr) = error1
     628           0 :       elist(last) = error2
     629           0 :       go to 20
     630           0 :    10 alist(maxerr) = a2
     631           0 :       alist(last) = a1
     632           0 :       blist(last) = b1
     633           0 :       rlist(maxerr) = area2
     634           0 :       rlist(last) = area1
     635           0 :       elist(maxerr) = error2
     636           0 :       elist(last) = error1
     637             :       !
     638             :       !           call subroutine dqpsrt to maintain the descending ordering
     639             :       !           in the list of error estimates and select the subinterval
     640             :       !           with the largest error estimate (to be bisected next).
     641             :       !
     642           0 :    20 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
     643             :       ! ***jump out of do-loop
     644           0 :       if(ier.ne.0.or.errsum.le.errbnd) go to 40
     645             :     end do
     646             :     !
     647             :     !           compute final result.
     648             :     !           ---------------------
     649             :     !
     650           0 :  40 result = 0.0_f
     651           0 :     do k=1,last
     652           0 :       result = result+rlist(k)
     653             :     end do
     654           0 :     abserr = errsum
     655           0 :  60 if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1)
     656           0 :     if(keyf.eq.1) neval = 30*neval+15
     657             : 999 return
     658             :   end subroutine dqage
     659             : 
     660             : 
     661             :   !!***begin prologue  dqagi
     662             :   !!***date written   800101   (yymmdd)
     663             :   !!***revision date  130319   (yymmdd)
     664             :   !!***category no.  h2a3a1,h2a4a1
     665             :   !!***keywords  automatic integrator, infinite intervals,
     666             :   !!             general-purpose, transformation, extrapolation,
     667             :   !!             globally adaptive
     668             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
     669             :   !!           de doncker,elise,appl. math. & progr. div. -k.u.leuven
     670             :   !!***purpose  the routine calculates an approximation result to a given
     671             :   !!            integral   i = integral of f over (bound,+infinity)
     672             :   !!            or i = integral of f over (-infinity,bound)
     673             :   !!            or i = integral of f over (-infinity,+infinity)
     674             :   !!            hopefully satisfying following claim for accuracy
     675             :   !!            abs(i-result).le.max(epsabs,epsrel*abs(i)).
     676             :   !!***description
     677             :   !!
     678             :   !!        integration over infinite intervals
     679             :   !!        standard fortran subroutine
     680             :   !!
     681             :   !!        parameters
     682             :   !!         on entry
     683             :   !!            fx     - double precision
     684             :   !!                     function subprogram defining the integrand
     685             :   !!                     function f(x). the actual name for f needs to be
     686             :   !!                     declared e x t e r n a l in the driver program.
     687             :   !!
     688             :   !!            fx_vars- structure containing variables need for integration
     689             :   !!                     specific to fractal meanfield scattering code
     690             :   !!
     691             :   !!            bound  - double precision
     692             :   !!                     finite bound of integration range
     693             :   !!                     (has no meaning if interval is doubly-infinite)
     694             :   !!
     695             :   !!            inf    - integer
     696             :   !!                     indicating the kind of integration range involved
     697             :   !!                     inf = 1 corresponds to  (bound,+infinity),
     698             :   !!                     inf = -1            to  (-infinity,bound),
     699             :   !!                     inf = 2             to (-infinity,+infinity).
     700             :   !!
     701             :   !!            epsabs - double precision
     702             :   !!                     absolute accuracy requested
     703             :   !!            epsrel - double precision
     704             :   !!                     relative accuracy requested
     705             :   !!                     if  epsabs.le.0
     706             :   !!                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
     707             :   !!                     the routine will end with ier = 6.
     708             :   !!
     709             :   !!
     710             :   !!         on return
     711             :   !!            result - double precision
     712             :   !!                     approximation to the integral
     713             :   !!
     714             :   !!            abserr - double precision
     715             :   !!                     estimate of the modulus of the absolute error,
     716             :   !!                     which should equal or exceed abs(i-result)
     717             :   !!
     718             :   !!            neval  - integer
     719             :   !!                     number of integrand evaluations
     720             :   !!
     721             :   !!            ier    - integer
     722             :   !!                     ier = 0 normal and reliable termination of the
     723             :   !!                             routine. it is assumed that the requested
     724             :   !!                             accuracy has been achieved.
     725             :   !!                   - ier.gt.0 abnormal termination of the routine. the
     726             :   !!                             estimates for result and error are less
     727             :   !!                             reliable. it is assumed that the requested
     728             :   !!                             accuracy has not been achieved.
     729             :   !!            error messages
     730             :   !!                     ier = 1 maximum number of subdivisions allowed
     731             :   !!                             has been achieved. one can allow more
     732             :   !!                             subdivisions by increasing the value of
     733             :   !!                             limit (and taking the according dimension
     734             :   !!                             adjustments into account). however, if
     735             :   !!                             this yields no improvement it is advised
     736             :   !!                             to analyze the integrand in order to
     737             :   !!                             determine the integration difficulties. if
     738             :   !!                             the position of a local difficulty can be
     739             :   !!                             determined (e.g. singularity,
     740             :   !!                             discontinuity within the interval) one
     741             :   !!                             will probably gain from splitting up the
     742             :   !!                             interval at this point and calling the
     743             :   !!                             integrator on the subranges. if possible,
     744             :   !!                             an appropriate special-purpose integrator
     745             :   !!                             should be used, which is designed for
     746             :   !!                             handling the type of difficulty involved.
     747             :   !!                         = 2 the occurrence of roundoff error is
     748             :   !!                             detected, which prevents the requested
     749             :   !!                             tolerance from being achieved.
     750             :   !!                             the error may be under-estimated.
     751             :   !!                         = 3 extremely bad integrand behaviour occurs
     752             :   !!                             at some points of the integration
     753             :   !!                             interval.
     754             :   !!                         = 4 the algorithm does not converge.
     755             :   !!                             roundoff error is detected in the
     756             :   !!                             extrapolation table.
     757             :   !!                             it is assumed that the requested tolerance
     758             :   !!                             cannot be achieved, and that the returned
     759             :   !!                             result is the best which can be obtained.
     760             :   !!                         = 5 the integral is probably divergent, or
     761             :   !!                             slowly convergent. it must be noted that
     762             :   !!                             divergence can occur with any other value
     763             :   !!                             of ier.
     764             :   !!                         = 6 the input is invalid, because
     765             :   !!                             (epsabs.le.0 and
     766             :   !!                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
     767             :   !!                              or limit.lt.1 or leniw.lt.limit*4.
     768             :   !!                             result, abserr, neval, last are set to
     769             :   !!                             zero. exept when limit or leniw is
     770             :   !!                             invalid, iwork(1), work(limit*2+1) and
     771             :   !!                             work(limit*3+1) are set to zero, work(1)
     772             :   !!                             is set to a and work(limit+1) to b.
     773             :   !!                         = 9 failure in sd1mach determining machine parameters
     774             :   !!
     775             :   !!          dimensioning parameters
     776             :   !!            limit - integer
     777             :   !!                    dimensioning parameter for iwork
     778             :   !!                    limit determines the maximum number of subintervals
     779             :   !!                    in the partition of the given integration interval
     780             :   !!                    (a,b), limit.ge.1.
     781             :   !!                    if limit.lt.1, the routine will end with ier = 6.
     782             :   !!
     783             :   !!            lenw  - integer
     784             :   !!                    dimensioning parameter for work
     785             :   !!                    lenw must be at least limit*4.
     786             :   !!                    if lenw.lt.limit*4, the routine will end
     787             :   !!                    with ier = 6.
     788             :   !!
     789             :   !!            last  - integer
     790             :   !!                    on return, last equals the number of subintervals
     791             :   !!                    produced in the subdivision process, which
     792             :   !!                    determines the number of significant elements
     793             :   !!                    actually in the work arrays.
     794             :   !!
     795             :   !!         work arrays
     796             :   !!            iwork - integer
     797             :   !!                    vector of dimension at least limit, the first
     798             :   !!                    k elements of which contain pointers
     799             :   !!                    to the error estimates over the subintervals,
     800             :   !!                    such that work(limit*3+iwork(1)),... ,
     801             :   !!                    work(limit*3+iwork(k)) form a decreasing
     802             :   !!                    sequence, with k = last if last.le.(limit/2+2), and
     803             :   !!                    k = limit+1-last otherwise
     804             :   !!
     805             :   !!            work  - double precision
     806             :   !!                    vector of dimension at least lenw
     807             :   !!                    on return
     808             :   !!                    work(1), ..., work(last) contain the left
     809             :   !!                     end points of the subintervals in the
     810             :   !!                     partition of (a,b),
     811             :   !!                    work(limit+1), ..., work(limit+last) contain
     812             :   !!                     the right end points,
     813             :   !!                    work(limit*2+1), ...,work(limit*2+last) contain the
     814             :   !!                     integral approximations over the subintervals,
     815             :   !!                    work(limit*3+1), ..., work(limit*3)
     816             :   !!                     contain the error estimates.
     817             :   !!***references  (none)
     818             :   !!***routines called  dqagie,xerror
     819             :   !!***end prologue  dqagi
     820             :   !!
     821           0 :   subroutine dqagi(fx,fx_vars,bound,inf,epsabs,epsrel,result,abserr,neval, &
     822           0 :                    ier,limit,lenw,last,iwork,work)
     823             : 
     824             :     ! Arguments
     825             :     interface
     826             :       function fx(centr, vars)
     827             :         use carma_precision_mod, only : f
     828             :         use adgaquad_types_mod
     829             :         real(kind=f), intent(in) :: centr
     830             :         type(adgaquad_vars_type), intent(inout) :: vars
     831             :         real(kind=f)             :: fx
     832             :       end function fx
     833             :     end interface
     834             :     type(adgaquad_vars_type) :: fx_vars
     835             :     real(kind=f) :: bound
     836             :     integer :: inf
     837             :     real(kind=f) :: epsabs 
     838             :     real(kind=f) :: epsrel
     839             :     real(kind=f) :: result
     840             :     real(kind=f) :: abserr
     841             :     integer :: neval
     842             :     integer :: ier
     843             :     integer :: limit  
     844             :     integer :: lenw  
     845             :     integer ::  last
     846             :     integer ::  iwork(limit)
     847             :     real(kind=f) :: work(lenw)
     848             : 
     849             :     ! Local declarations
     850             :     integer lvl,l1,l2,l3
     851             : 
     852             :     !
     853             :     !         check validity of limit and lenw.
     854             :     !
     855             :     !***first executable statement  dqagi
     856           0 :     ier = 6
     857           0 :     neval = 0
     858           0 :     last = 0
     859           0 :     result = 0.0_f
     860           0 :     abserr = 0.0_f
     861           0 :     if(limit.lt.1.or.lenw.lt.limit*4) go to 10
     862             :     !
     863             :     !         prepare call for dqagie.
     864             :     !
     865           0 :     l1 = limit+1
     866           0 :     l2 = limit+l1
     867           0 :     l3 = limit+l2
     868             : 
     869             :     call dqagie(fx,fx_vars,bound,inf,epsabs,epsrel,limit,result,abserr, &
     870           0 :                 neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last)
     871             :     !
     872             :     !         call error handler if necessary.
     873             :     !
     874           0 :     lvl = 0
     875           0 : 10  if(ier.eq.6) lvl = 1
     876           0 :     if(ier.ne.0) then
     877           0 :       write(*,*) "ERROR: abnormal return from dqagi"
     878           0 :       write(*,*) "       ifail=",ier,"  level=",lvl
     879             :     endif
     880           0 :     return
     881             :   end subroutine dqagi
     882             : 
     883             : 
     884             :   !!***begin prologue  dqagie
     885             :   !!***date written   800101   (yymmdd)
     886             :   !!***revision date  130319   (yymmdd)
     887             :   !!***category no.  h2a3a1,h2a4a1
     888             :   !!***keywords  automatic integrator, infinite intervals,
     889             :   !!             general-purpose, transformation, extrapolation,
     890             :   !!             globally adaptive
     891             :   !!***author  piessens,robert,appl. math & progr. div - k.u.leuven
     892             :   !!           de doncker,elise,appl. math & progr. div - k.u.leuven
     893             :   !!***purpose  the routine calculates an approximation result to a given
     894             :   !!            integral   i = integral of f over (bound,+infinity)
     895             :   !!            or i = integral of f over (-infinity,bound)
     896             :   !!            or i = integral of f over (-infinity,+infinity),
     897             :   !!            hopefully satisfying following claim for accuracy
     898             :   !!            abs(i-result).le.max(epsabs,epsrel*abs(i))
     899             :   !!***description
     900             :   !!
     901             :   !! integration over infinite intervals
     902             :   !! standard fortran subroutine
     903             :   !!
     904             :   !!            fx     - double precision
     905             :   !!                     function subprogram defining the integrand
     906             :   !!                     function f(x). the actual name for f needs to be
     907             :   !!                     declared e x t e r n a l in the driver program.
     908             :   !!
     909             :   !!            fx_vars- structure containing variables need for integration
     910             :   !!                     specific to fractal meanfield scattering code
     911             :   !!
     912             :   !!            bound  - double precision
     913             :   !!                     finite bound of integration range
     914             :   !!                     (has no meaning if interval is doubly-infinite)
     915             :   !!
     916             :   !!            inf    - double precision
     917             :   !!                     indicating the kind of integration range involved
     918             :   !!                     inf = 1 corresponds to  (bound,+infinity),
     919             :   !!                     inf = -1            to  (-infinity,bound),
     920             :   !!                     inf = 2             to (-infinity,+infinity).
     921             :   !!
     922             :   !!            epsabs - double precision
     923             :   !!                     absolute accuracy requested
     924             :   !!            epsrel - double precision
     925             :   !!                     relative accuracy requested
     926             :   !!                     if  epsabs.le.0
     927             :   !!                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
     928             :   !!                     the routine will end with ier = 6.
     929             :   !!
     930             :   !!            limit  - integer
     931             :   !!                     gives an upper bound on the number of subintervals
     932             :   !!                     in the partition of (a,b), limit.ge.1
     933             :   !!
     934             :   !!         on return
     935             :   !!           result - double precision
     936             :   !!                     approximation to the integral
     937             :   !!
     938             :   !!            abserr - double precision
     939             :   !!                     estimate of the modulus of the absolute error,
     940             :   !!                     which should equal or exceed abs(i-result)
     941             :   !!
     942             :   !!            neval  - integer
     943             :   !!                     number of integrand evaluations
     944             :   !!
     945             :   !!            ier    - integer
     946             :   !!                     ier = 0 normal and reliable termination of the
     947             :   !!                             routine. it is assumed that the requested
     948             :   !!                             accuracy has been achieved.
     949             :   !!                   - ier.gt.0 abnormal termination of the routine. the
     950             :   !!                             estimates for result and error are less
     951             :   !!                             reliable. it is assumed that the requested
     952             :   !!                             accuracy has not been achieved.
     953             :   !!            error messages
     954             :   !!                     ier = 1 maximum number of subdivisions allowed
     955             :   !!                             has been achieved. one can allow more
     956             :   !!                             subdivisions by increasing the value of
     957             :   !!                             limit (and taking the according dimension
     958             :   !!                             adjustments into account). however,if
     959             :   !!                             this yields no improvement it is advised
     960             :   !!                             to analyze the integrand in order to
     961             :   !!                             determine the integration difficulties.
     962             :   !!                             if the position of a local difficulty can
     963             :   !!                             be determined (e.g. singularity,
     964             :   !!                             discontinuity within the interval) one
     965             :   !!                             will probably gain from splitting up the
     966             :   !!                             interval at this point and calling the
     967             :   !!                             integrator on the subranges. if possible,
     968             :   !!                             an appropriate special-purpose integrator
     969             :   !!                             should be used, which is designed for
     970             :   !!                             handling the type of difficulty involved.
     971             :   !!                         = 2 the occurrence of roundoff error is
     972             :   !!                             detected, which prevents the requested
     973             :   !!                             tolerance from being achieved.
     974             :   !!                             the error may be under-estimated.
     975             :   !!                         = 3 extremely bad integrand behaviour occurs
     976             :   !!                             at some points of the integration
     977             :   !!                             interval.
     978             :   !!                         = 4 the algorithm does not converge.
     979             :   !!                             roundoff error is detected in the
     980             :   !!                             extrapolation table.
     981             :   !!                             it is assumed that the requested tolerance
     982             :   !!                             cannot be achieved, and that the returned
     983             :   !!                             result is the best which can be obtained.
     984             :   !!                         = 5 the integral is probably divergent, or
     985             :   !!                             slowly convergent. it must be noted that
     986             :   !!                             divergence can occur with any other value
     987             :   !!                             of ier.
     988             :   !!                         = 6 the input is invalid, because
     989             :   !!                             (epsabs.le.0 and
     990             :   !!                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
     991             :   !!                             result, abserr, neval, last, rlist(1),
     992             :   !!                             elist(1) and iord(1) are set to zero.
     993             :   !!                             alist(1) and blist(1) are set to 0
     994             :   !!                             and 1 respectively.
     995             :   !!                         = 9 failure in sd1mach determining machine parameters
     996             :   !!
     997             :   !!            alist  - double precision
     998             :   !!                     vector of dimension at least limit, the first
     999             :   !!                      last  elements of which are the left
    1000             :   !!                     end points of the subintervals in the partition
    1001             :   !!                     of the transformed integration range (0,1).
    1002             :   !!
    1003             :   !!            blist  - double precision
    1004             :   !!                     vector of dimension at least limit, the first
    1005             :   !!                      last  elements of which are the right
    1006             :   !!                     end points of the subintervals in the partition
    1007             :   !!                     of the transformed integration range (0,1).
    1008             :   !!
    1009             :   !!            rlist  - double precision
    1010             :   !!                     vector of dimension at least limit, the first
    1011             :   !!                      last  elements of which are the integral
    1012             :   !!                     approximations on the subintervals
    1013             :   !!
    1014             :   !!            elist  - double precision
    1015             :   !!                     vector of dimension at least limit,  the first
    1016             :   !!                     last elements of which are the moduli of the
    1017             :   !!                     absolute error estimates on the subintervals
    1018             :   !!
    1019             :   !!            iord   - integer
    1020             :   !!                     vector of dimension limit, the first k
    1021             :   !!                     elements of which are pointers to the
    1022             :   !!                     error estimates over the subintervals,
    1023             :   !!                     such that elist(iord(1)), ..., elist(iord(k))
    1024             :   !!                     form a decreasing sequence, with k = last
    1025             :   !!                     if last.le.(limit/2+2), and k = limit+1-last
    1026             :   !!                     otherwise
    1027             :   !!
    1028             :   !!            last   - integer
    1029             :   !!                     number of subintervals actually produced
    1030             :   !!                     in the subdivision process
    1031             :   !!
    1032             :   !!***references  (none)
    1033             :   !!***routines called  sd1mach,dqelg,dqk15i,dqpsrt
    1034             :   !!***end prologue  dqagie
    1035           0 :   subroutine dqagie(fx,fx_vars,bound,inf,epsabs,epsrel,limit,result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)
    1036             : 
    1037             :     ! Arguments
    1038             :     interface
    1039             :       function fx(centr, vars)
    1040             :         use carma_precision_mod, only : f
    1041             :         use adgaquad_types_mod
    1042             :         real(kind=f), intent(in) :: centr
    1043             :         type(adgaquad_vars_type), intent(inout) :: vars
    1044             :         real(kind=f)             :: fx
    1045             :       end function fx
    1046             :     end interface
    1047             :     type(adgaquad_vars_type) :: fx_vars
    1048             :     real(kind=f) :: bound
    1049             :     integer :: inf
    1050             :     real(kind=f) :: epsabs
    1051             :     real(kind=f) :: epsrel
    1052             :     integer :: limit
    1053             :     real(kind=f) :: result
    1054             :     real(kind=f) :: abserr
    1055             :     integer :: neval
    1056             :     integer :: ier
    1057             :     real(kind=f) :: alist(limit)
    1058             :     real(kind=f) :: blist(limit)
    1059             :     real(kind=f) :: rlist(limit)
    1060             :     real(kind=f) :: elist(limit)
    1061             :     integer :: iord(limit)
    1062             :     integer :: last
    1063             :     
    1064             :     ! Local declartions
    1065             :     real(kind=f) :: abseps, area, area1, area12, area2
    1066             :     real(kind=f) :: a1, a2, b1, b2, correc
    1067             :     real(kind=f) :: defabs, defab1, defab2
    1068             :     real(kind=f) :: dmax1, dres, epmach, erlarg, erlast
    1069             :     real(kind=f) :: errbnd, errmax, error1, error2, erro12, errsum
    1070             :     real(kind=f) :: ertest, oflow, resabs, reseps, res3la(3), rlist2(52)
    1071             :     real(kind=f) :: small, uflow, boun
    1072             :       
    1073             :     integer :: id, ierro, iroff1, iroff2, iroff3, jupbnd, k, ksgn
    1074             :     integer :: ktmin, maxerr, nres, nrmax, numrl2
    1075             :     logical :: extrap, noext
    1076             : 
    1077             : 
    1078             :     !
    1079             :     !            the dimension of rlist2 is determined by the value of
    1080             :     !            limexp in subroutine dqelg.
    1081             :     !
    1082             :     !            list of major variables
    1083             :     !            -----------------------
    1084             :     !
    1085             :     !           alist     - list of left end points of all subintervals
    1086             :     !                       considered up to now
    1087             :     !           blist     - list of right end points of all subintervals
    1088             :     !                       considered up to now
    1089             :     !           rlist(i)  - approximation to the integral over
    1090             :     !                       (alist(i),blist(i))
    1091             :     !           rlist2    - array of dimension at least (limexp+2),
    1092             :     !                       containing the part of the epsilon table
    1093             :     !                       wich is still needed for further computations
    1094             :     !           elist(i)  - error estimate applying to rlist(i)
    1095             :     !           maxerr    - pointer to the interval with largest error
    1096             :     !                       estimate
    1097             :     !           errmax    - elist(maxerr)
    1098             :     !           erlast    - error on the interval currently subdivided
    1099             :     !                       (before that subdivision has taken place)
    1100             :     !           area      - sum of the integrals over the subintervals
    1101             :     !           errsum    - sum of the errors over the subintervals
    1102             :     !           errbnd    - requested accuracy max(epsabs,epsrel*
    1103             :     !                       abs(result))
    1104             :     !           *****1    - variable for the left subinterval
    1105             :     !           *****2    - variable for the right subinterval
    1106             :     !           last      - index for subdivision
    1107             :     !           nres      - number of calls to the extrapolation routine
    1108             :     !           numrl2    - number of elements currently in rlist2. if an
    1109             :     !                       appropriate approximation to the compounded
    1110             :     !                       integral has been obtained, it is put in
    1111             :     !                       rlist2(numrl2) after numrl2 has been increased
    1112             :     !                       by one.
    1113             :     !           small     - length of the smallest interval considered up
    1114             :     !                       to now, multiplied by 1.5
    1115             :     !           erlarg    - sum of the errors over the intervals larger
    1116             :     !                       than the smallest interval considered up to now
    1117             :     !           extrap    - logical variable denoting that the routine
    1118             :     !                       is attempting to perform extrapolation. i.e.
    1119             :     !                       before subdividing the smallest interval we
    1120             :     !                       try to decrease the value of erlarg.
    1121             :     !           noext     - logical variable denoting that extrapolation
    1122             :     !                       is no longer allowed (true-value)
    1123             :     !
    1124             :     !            machine dependent constants
    1125             :     !            ---------------------------
    1126             :     !
    1127             :     !           epmach is the largest relative spacing.
    1128             :     !           uflow is the smallest positive magnitude.
    1129             :     !           oflow is the largest positive magnitude.
    1130             :     !
    1131             :     !***first executable statement  dqagie
    1132             : 
    1133           0 :     call sd1mach(4,epmach,ier)
    1134           0 :     if(ier.eq.9) return
    1135             :     
    1136             :     !epmach = d1mach(4)
    1137             :     !
    1138             :     !           test on validity of parameters
    1139             :     !           -----------------------------
    1140             :     !
    1141           0 :     ier = 0
    1142           0 :     neval = 0
    1143           0 :     last = 0
    1144           0 :     result = 0.0_f
    1145           0 :     abserr = 0.0_f
    1146           0 :     alist(1) = 0.0_f
    1147           0 :     blist(1) = 0.1e1_f
    1148           0 :     rlist(1) = 0.0_f
    1149           0 :     elist(1) = 0.0_f
    1150           0 :     iord(1) = 0
    1151           0 :     if(epsabs.le.0.0_f.and.epsrel.lt.dmax1(0.5e2_f*epmach,0.5e-28_f)) ier = 6
    1152           0 :     if(ier.eq.6) go to 999
    1153             :     !
    1154             :     !
    1155             :     !           first approximation to the integral
    1156             :     !           -----------------------------------
    1157             :     !
    1158             :     !           determine the interval to be mapped onto (0,1).
    1159             :     !           if inf = 2 the integral is computed as i = i1+i2, where
    1160             :     !           i1 = integral of f over (-infinity,0),
    1161             :     !           i2 = integral of f over (0,+infinity).
    1162             :     !
    1163           0 :     boun = bound
    1164           0 :     if(inf.eq.2) boun = 0.0_f
    1165           0 :     call dqk15i(fx,fx_vars,boun,inf,0.0_f,0.1e1_f,result,abserr,defabs,resabs,ier)
    1166           0 :     if(ier.eq.9) return
    1167             :     !
    1168             :     !           test on accuracy
    1169             :     !
    1170           0 :     last = 1
    1171           0 :     rlist(1) = result
    1172           0 :     elist(1) = abserr
    1173           0 :     iord(1) = 1
    1174           0 :     dres = dabs(result)
    1175           0 :     errbnd = dmax1(epsabs,epsrel*dres)
    1176           0 :     if(abserr.le.1.0e2_f*epmach*defabs.and.abserr.gt.errbnd) ier = 2
    1177           0 :     if(limit.eq.1) ier = 1
    1178           0 :     if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or.abserr.eq.0.0_f) go to 130
    1179             :     !
    1180             :     !           initialization
    1181             :     !           --------------
    1182             :     !
    1183           0 :     call sd1mach(1,uflow,ier)
    1184           0 :     if(ier.eq.9) return
    1185           0 :     call sd1mach(2,oflow,ier)
    1186           0 :     if(ier.eq.9) return
    1187             : 
    1188             :     !uflow = d1mach(1)
    1189             :     !oflow = d1mach(2)
    1190           0 :     rlist2(1) = result
    1191           0 :     errmax = abserr
    1192           0 :     maxerr = 1
    1193           0 :     area = result
    1194           0 :     errsum = abserr
    1195           0 :     abserr = oflow
    1196           0 :     nrmax = 1
    1197           0 :     nres = 0
    1198           0 :     ktmin = 0
    1199           0 :     numrl2 = 2
    1200           0 :     extrap = .false.
    1201           0 :     noext = .false.
    1202           0 :     ierro = 0
    1203           0 :     iroff1 = 0
    1204           0 :     iroff2 = 0
    1205           0 :     iroff3 = 0
    1206           0 :     ksgn = -1
    1207           0 :     if(dres.ge.(0.1e1_f-0.5e2_f*epmach)*defabs) ksgn = 1
    1208             :     !
    1209             :     !           main do-loop
    1210             :     !           ------------
    1211             :     !
    1212           0 :     do 90 last = 2,limit
    1213             :       !
    1214             :       !           bisect the subinterval with nrmax-th largest error estimate.
    1215             :       !
    1216           0 :       a1 = alist(maxerr)
    1217           0 :       b1 = 0.5_f*(alist(maxerr)+blist(maxerr))
    1218           0 :       a2 = b1
    1219           0 :       b2 = blist(maxerr)
    1220           0 :       erlast = errmax
    1221           0 :       call dqk15i(fx,fx_vars,boun,inf,a1,b1,area1,error1,resabs,defab1,ier)
    1222           0 :       if(ier.eq.9) return
    1223           0 :       call dqk15i(fx,fx_vars,boun,inf,a2,b2,area2,error2,resabs,defab2,ier)
    1224           0 :       if(ier.eq.9) return
    1225             :       !
    1226             :       !           improve previous approximations to integral
    1227             :       !           and error and test for accuracy.
    1228             :       !
    1229           0 :       area12 = area1+area2
    1230           0 :       erro12 = error1+error2
    1231           0 :       errsum = errsum+erro12-errmax
    1232           0 :       area = area+area12-rlist(maxerr)
    1233           0 :       if(defab1.eq.error1.or.defab2.eq.error2)go to 15
    1234           0 :       if(dabs(rlist(maxerr)-area12).gt.0.1e-4_f*dabs(area12).or.erro12.lt.0.99_f*errmax) go to 10
    1235           0 :       if(extrap) iroff2 = iroff2+1
    1236           0 :       if(.not.extrap) iroff1 = iroff1+1
    1237           0 :  10   if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
    1238           0 :  15   rlist(maxerr) = area1
    1239           0 :       rlist(last) = area2
    1240           0 :       errbnd = dmax1(epsabs,epsrel*dabs(area))
    1241             :       !
    1242             :       !           test for roundoff error and eventually set error flag.
    1243             :       !
    1244           0 :       if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
    1245           0 :       if(iroff2.ge.5) ierro = 3
    1246             :       !
    1247             :       !           set error flag in the case that the number of
    1248             :       !           subintervals equals limit.
    1249             :       !
    1250           0 :       if(last.eq.limit) ier = 1
    1251             :       !
    1252             :       !           set error flag in the case of bad integrand behaviour
    1253             :       !           at some points of the integration range.
    1254             :       !
    1255           0 :       if(dmax1(dabs(a1),dabs(b2)).le.(0.1e1_f+0.1e3_f*epmach)*(dabs(a2)+0.1e4_f*uflow)) ier = 4
    1256             :       !
    1257             :       !           append the newly-created intervals to the list.
    1258             :       !
    1259           0 :       if(error2.gt.error1) go to 20
    1260           0 :       alist(last) = a2
    1261           0 :       blist(maxerr) = b1
    1262           0 :       blist(last) = b2
    1263           0 :       elist(maxerr) = error1
    1264           0 :       elist(last) = error2
    1265           0 :       go to 30
    1266           0 :  20   alist(maxerr) = a2
    1267           0 :       alist(last) = a1
    1268           0 :       blist(last) = b1
    1269           0 :       rlist(maxerr) = area2
    1270           0 :       rlist(last) = area1
    1271           0 :       elist(maxerr) = error2
    1272           0 :       elist(last) = error1
    1273             :       !
    1274             :       !           call subroutine dqpsrt to maintain the descending ordering
    1275             :       !           in the list of error estimates and select the subinterval
    1276             :       !           with nrmax-th largest error estimate (to be bisected next).
    1277             :       !
    1278           0 :  30   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
    1279           0 :       if(errsum.le.errbnd) go to 115
    1280           0 :       if(ier.ne.0) go to 100
    1281           0 :       if(last.eq.2) go to 80
    1282           0 :       if(noext) go to 90
    1283           0 :       erlarg = erlarg-erlast
    1284           0 :       if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12
    1285           0 :       if(extrap) go to 40
    1286             :       !
    1287             :       !           test whether the interval to be bisected next is the
    1288             :       !           smallest interval.
    1289             :       !
    1290           0 :       if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
    1291           0 :       extrap = .true.
    1292           0 :       nrmax = 2
    1293           0 :  40   if(ierro.eq.3.or.erlarg.le.ertest) go to 60
    1294             :       !
    1295             :       !           the smallest interval has the largest error.
    1296             :       !           before bisecting decrease the sum of the errors over the
    1297             :       !           larger intervals (erlarg) and perform extrapolation.
    1298             :       !
    1299           0 :       id = nrmax
    1300           0 :       jupbnd = last
    1301           0 :       if(last.gt.(2+limit/2)) jupbnd = limit+3-last
    1302           0 :       do k = id,jupbnd
    1303           0 :         maxerr = iord(nrmax)
    1304           0 :         errmax = elist(maxerr)
    1305           0 :         if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90
    1306           0 :         nrmax = nrmax+1
    1307             :       end do
    1308             :       !
    1309             :       !           perform extrapolation.
    1310             :       !
    1311           0 :  60   numrl2 = numrl2+1
    1312           0 :       rlist2(numrl2) = area
    1313           0 :       call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres, ier)
    1314           0 :       if(ier.eq.9) return
    1315           0 :       ktmin = ktmin+1
    1316           0 :       if(ktmin.gt.5.and.abserr.lt.0.1e-2_f*errsum) ier = 5
    1317           0 :       if(abseps.ge.abserr) go to 70
    1318           0 :       ktmin = 0
    1319           0 :       abserr = abseps
    1320           0 :       result = reseps
    1321           0 :       correc = erlarg
    1322           0 :       ertest = dmax1(epsabs,epsrel*dabs(reseps))
    1323           0 :       if(abserr.le.ertest) go to 100
    1324             :       !
    1325             :       !            prepare bisection of the smallest interval.
    1326             :       !
    1327           0 :  70   if(numrl2.eq.1) noext = .true.
    1328           0 :       if(ier.eq.5) go to 100
    1329           0 :       maxerr = iord(1)
    1330           0 :       errmax = elist(maxerr)
    1331           0 :       nrmax = 1
    1332           0 :       extrap = .false.
    1333           0 :       small = small*0.5_f
    1334           0 :       erlarg = errsum
    1335           0 :       go to 90
    1336           0 :  80   small = 0.375_f
    1337           0 :       erlarg = errsum
    1338           0 :       ertest = errbnd
    1339           0 :       rlist2(2) = area
    1340           0 :  90 continue
    1341             :     !
    1342             :     !           set final result and error estimate.
    1343             :     !           ------------------------------------
    1344             :     !
    1345           0 : 100 if(abserr.eq.oflow) go to 115
    1346           0 :     if((ier+ierro).eq.0) go to 110
    1347           0 :     if(ierro.eq.3) abserr = abserr+correc
    1348           0 :     if(ier.eq.0) ier = 3
    1349           0 :     if(result.ne.0.0_f.and.area.ne.0.0_f)go to 105
    1350           0 :     if(abserr.gt.errsum)go to 115
    1351           0 :     if(area.eq.0.0_f) go to 130
    1352           0 :     go to 110
    1353           0 : 105 if(abserr/dabs(result).gt.errsum/dabs(area))go to 115
    1354             :     !
    1355             :     !           test on divergence
    1356             :     !
    1357           0 : 110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le.defabs*0.1e-1_f) go to 130
    1358           0 :     if(0.1e-1_f.gt.(result/area).or.(result/area).gt.0.1e3_f.or.errsum.gt.dabs(area)) ier = 6
    1359           0 :     go to 130
    1360             :     !
    1361             :     !           compute global integral sum.
    1362             :     !
    1363           0 : 115 result = 0.0_f
    1364           0 :     do k = 1,last
    1365           0 :       result = result+rlist(k)
    1366             :     end do
    1367           0 :     abserr = errsum
    1368           0 : 130 neval = 30*last-15
    1369           0 :     if(inf.eq.2) neval = 2*neval
    1370           0 :     if(ier.gt.2) ier=ier-1
    1371             : 999 return
    1372             :   end subroutine dqagie
    1373             : 
    1374             : 
    1375             :   !!***begin prologue  dqk15
    1376             :   !!***date written   800101   (yymmdd)
    1377             :   !!***revision date  130319   (yymmdd)
    1378             :   !!***category no.  h2a1a2
    1379             :   !!***keywords  15-point gauss-kronrod rules
    1380             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    1381             :   !!           de doncker,elise,appl. math. & progr. div - k.u.leuven
    1382             :   !!***purpose  to compute i = integral of f over (a,b), with error
    1383             :   !!                           estimate
    1384             :   !!                       j = integral of abs(f) over (a,b)
    1385             :   !!***description
    1386             :   !!
    1387             :   !!           integration rules
    1388             :   !!           standard fortran subroutine
    1389             :   !!           double precision version
    1390             :   !!
    1391             :   !!           parameters
    1392             :   !!            on entry
    1393             :   !!              fx     - double precision
    1394             :   !!                       function subprogram defining the integrand
    1395             :   !!                       function f(x). the actual name for f needs to be
    1396             :   !!                       declared e x t e r n a l in the calling program.
    1397             :   !!
    1398             :   !!              fx_vars- structure containing variables need for integration
    1399             :   !!                     specific to fractal meanfield scattering code!
    1400             :   !!
    1401             :   !!              a      - double precision
    1402             :   !!                       lower limit of integration
    1403             :   !!
    1404             :   !!              b      - double precision
    1405             :   !!                       upper limit of integration
    1406             :   !!
    1407             :   !!            on return
    1408             :   !!              result - double precision
    1409             :   !!                       approximation to the integral i
    1410             :   !!                       result is computed by applying the 15-point
    1411             :   !!                       kronrod rule (resk) obtained by optimal addition
    1412             :   !!                       of abscissae to the7-point gauss rule(resg).
    1413             :   !!
    1414             :   !!              abserr - double precision
    1415             :   !!                       estimate of the modulus of the absolute error,
    1416             :   !!                       which should not exceed abs(i-result)
    1417             :   !!
    1418             :   !!              resabs - double precision
    1419             :   !!                       approximation to the integral j
    1420             :   !!
    1421             :   !!              resasc - double precision
    1422             :   !!                       approximation to the integral of abs(f-i/(b-a))
    1423             :   !!                       over (a,b)
    1424             :   !!
    1425             :   !!              ier    - integer
    1426             :   !!                       ier = 0 normal and reliable termination of the
    1427             :   !!                               routine. it is assumed that the requested
    1428             :   !!                               accuracy has been achieved.
    1429             :   !!                       ier.gt.0 abnormal termination of the routine. the
    1430             :   !!                                estimates for result and error are less
    1431             :   !!                                reliable. it is assumed that the requested
    1432             :   !!                                accuracy has not been achieved.
    1433             :   !!
    1434             :   !!***references  (none)
    1435             :   !!***routines called  sd1mach
    1436             :   !!***end prologue  dqk15  
    1437             :   !!
    1438           0 :   subroutine dqk15(fx,fx_vars,a,b,result,abserr,resabs,resasc,ier)
    1439             : 
    1440             :     ! Arguments
    1441             :     interface
    1442             :       function fx(centr, vars)
    1443             :         use carma_precision_mod, only : f
    1444             :         use adgaquad_types_mod
    1445             :         real(kind=f), intent(in) :: centr
    1446             :         type(adgaquad_vars_type), intent(inout) :: vars
    1447             :         real(kind=f)             :: fx
    1448             :       end function fx
    1449             :     end interface
    1450             :     type(adgaquad_vars_type) :: fx_vars
    1451             :     real(kind=f) :: a
    1452             :     real(kind=f) :: b 
    1453             :     real(kind=f) :: result
    1454             :     real(kind=f) :: abserr
    1455             :     real(kind=f) :: resabs
    1456             :     real(kind=f) :: resasc
    1457             :     integer :: ier
    1458             : 
    1459             :     ! Local Declarations
    1460             :     real(kind=f) :: absc, centr, dabs, dhlgth, dmax2, dmin1
    1461             :     real(kind=f) :: epmach, fc, fsum, fval1, fval2, fv1(7), fv2(7), hlgth
    1462             :     real(kind=f) :: resg, resk, reskh, uflow, wg(4), wgk(8), xgk(8)
    1463             :     integer :: j,jtw,jtwm1
    1464             : 
    1465             :     !
    1466             :     !
    1467             :     !           the abscissae and weights are given for the interval (-1,1).
    1468             :     !           because of symmetry only the positive abscissae and their
    1469             :     !           corresponding weights are given.
    1470             :     !
    1471             :     !           xgk    - abscissae of the 15-point kronrod rule
    1472             :     !                    xgk(2), xgk(4), ...  abscissae of the 7-point
    1473             :     !                    gauss rule
    1474             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    1475             :     !                    added to the 7-point gauss rule
    1476             :     !
    1477             :     !           wgk    - weights of the 15-point kronrod rule
    1478             :     !
    1479             :     !           wg     - weights of the 7-point gauss rule
    1480             :     !
    1481             :     !
    1482             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    1483             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    1484             :     ! bell labs, nov. 1981.
    1485             :     !
    1486             :     data wg  (  1) / 0.129484966168869693270611432679082_f /
    1487             :     data wg  (  2) / 0.279705391489276667901467771423780_f /
    1488             :     data wg  (  3) / 0.381830050505118944950369775488975_f /
    1489             :     data wg  (  4) / 0.417959183673469387755102040816327_f /
    1490             : 
    1491             :     data xgk (  1) / 0.991455371120812639206854697526329_f /
    1492             :     data xgk (  2) / 0.949107912342758524526189684047851_f /
    1493             :     data xgk (  3) / 0.864864423359769072789712788640926_f /
    1494             :     data xgk (  4) / 0.741531185599394439863864773280788_f /
    1495             :     data xgk (  5) / 0.586087235467691130294144838258730_f /
    1496             :     data xgk (  6) / 0.405845151377397166906606412076961_f /
    1497             :     data xgk (  7) / 0.207784955007898467600689403773245_f /
    1498             :     data xgk (  8) / 0.000000000000000000000000000000000_f /
    1499             : 
    1500             :     data wgk (  1) / 0.022935322010529224963732008058970_f /
    1501             :     data wgk (  2) / 0.063092092629978553290700663189204_f /
    1502             :     data wgk (  3) / 0.104790010322250183839876322541518_f /
    1503             :     data wgk (  4) / 0.140653259715525918745189590510238_f /
    1504             :     data wgk (  5) / 0.169004726639267902826583426598550_f /
    1505             :     data wgk (  6) / 0.190350578064785409913256402421014_f /
    1506             :     data wgk (  7) / 0.204432940075298892414161999234649_f /
    1507             :     data wgk (  8) / 0.209482141084727828012999174891714_f /
    1508             : 
    1509             :     !
    1510             :     !
    1511             :     !           list of major variables
    1512             :     !           -----------------------
    1513             :     !
    1514             :     !           centr  - mid point of the interval
    1515             :     !           hlgth  - half-length of the interval
    1516             :     !           absc   - abscissa
    1517             :     !           fval*  - function value
    1518             :     !           resg   - result of the 7-point gauss formula
    1519             :     !           resk   - result of the 15-point kronrod formula
    1520             :     !           reskh  - approximation to the mean value of f over (a,b),
    1521             :     !                    i.e. to i/(b-a)
    1522             :     !
    1523             :     !           machine dependent constants
    1524             :     !           ---------------------------
    1525             :     !
    1526             :     !           epmach is the largest relative spacing.
    1527             :     !           uflow is the smallest positive magnitude.
    1528             :     !
    1529             :     !***first executable statement  dqk15
    1530             :     !epmach = d1mach(4)
    1531             :     !uflow = d1mach(1)
    1532             : 
    1533           0 :     call sd1mach(4,epmach,ier)
    1534           0 :     if(ier.eq.9) return
    1535           0 :     call sd1mach(1,uflow,ier)
    1536           0 :     if(ier.eq.9) return
    1537             : 
    1538           0 :     centr = 0.5_f*(a+b)
    1539           0 :     hlgth = 0.5_f*(b-a)
    1540           0 :     dhlgth = dabs(hlgth)
    1541             :     !
    1542             :     !           compute the 15-point kronrod approximation to
    1543             :     !           the integral, and estimate the absolute error.
    1544             :     !
    1545           0 :     fc = fx(centr,fx_vars)
    1546           0 :     resg = fc*wg(4)
    1547           0 :     resk = fc*wgk(8)
    1548           0 :     resabs = dabs(resk)
    1549           0 :     do j=1,3
    1550           0 :       jtw = j*2
    1551           0 :       absc = hlgth*xgk(jtw)
    1552           0 :       fval1 = fx(centr-absc,fx_vars)
    1553           0 :       fval2 = fx(centr+absc,fx_vars)
    1554           0 :       fv1(jtw) = fval1
    1555           0 :       fv2(jtw) = fval2
    1556           0 :       fsum = fval1+fval2
    1557           0 :       resg = resg+wg(j)*fsum
    1558           0 :       resk = resk+wgk(jtw)*fsum
    1559           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    1560             :     end do
    1561           0 :     do j = 1,4
    1562           0 :       jtwm1 = j*2-1
    1563           0 :       absc = hlgth*xgk(jtwm1)
    1564           0 :       fval1 = fx(centr-absc,fx_vars)
    1565           0 :       fval2 = fx(centr+absc,fx_vars)
    1566           0 :       fv1(jtwm1) = fval1
    1567           0 :       fv2(jtwm1) = fval2
    1568           0 :       fsum = fval1+fval2
    1569           0 :       resk = resk+wgk(jtwm1)*fsum
    1570           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    1571             :     end do
    1572           0 :     reskh = resk*0.5_f
    1573           0 :     resasc = wgk(8)*dabs(fc-reskh)
    1574           0 :     do j=1,7
    1575           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    1576             :     end do
    1577           0 :     result = resk*hlgth
    1578           0 :     resabs = resabs*dhlgth
    1579           0 :     resasc = resasc*dhlgth
    1580           0 :     abserr = dabs((resk-resg)*hlgth)
    1581           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0.0_f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    1582           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    1583             : 999 return
    1584             :   end subroutine dqk15
    1585             : 
    1586             :   !!
    1587             :   !!***begin prologue  dqk21
    1588             :   !!***date written   800101   (yymmdd)
    1589             :   !!***revision date  130319   (yymmdd)
    1590             :   !!***category no.  h2a1a2
    1591             :   !!***keywords  21-point gauss-kronrod rules
    1592             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    1593             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    1594             :   !!***purpose  to compute i = integral of f over (a,b), with error
    1595             :   !!                           estimate
    1596             :   !!                       j = integral of abs(f) over (a,b)
    1597             :   !!***description
    1598             :   !!
    1599             :   !!           integration rules
    1600             :   !!           standard fortran subroutine
    1601             :   !!           double precision version
    1602             :   !!
    1603             :   !!           parameters
    1604             :   !!            on entry
    1605             :   !!              fx     - double precision
    1606             :   !!                       function subprogram defining the integrand
    1607             :   !!                       function f(x). the actual name for f needs to be
    1608             :   !!                       declared e x t e r n a l in the driver program.
    1609             :   !!
    1610             :   !!            fx_vars- structure containing variables need for integration
    1611             :   !!                     specific to fractal meanfield scattering code
    1612             :   !!              a      - double precision
    1613             :   !!                       lower limit of integration
    1614             :   !!
    1615             :   !!              b      - double precision
    1616             :   !!                       upper limit of integration
    1617             :   !!
    1618             :   !!            on return
    1619             :   !!              result - double precision
    1620             :   !!                       approximation to the integral i
    1621             :   !!                       result is computed by applying the 21-point
    1622             :   !!                       kronrod rule (resk) obtained by optimal addition
    1623             :   !!                       of abscissae to the 10-point gauss rule (resg).
    1624             :   !!
    1625             :   !!              abserr - double precision
    1626             :   !!                       estimate of the modulus of the absolute error,
    1627             :   !!                       which should not exceed abs(i-result)
    1628             :   !!
    1629             :   !!              resabs - double precision
    1630             :   !!                       approximation to the integral j
    1631             :   !!
    1632             :   !!              resasc - double precision
    1633             :   !!                       approximation to the integral of abs(f-i/(b-a))
    1634             :   !!                       over (a,b)
    1635             :   !!
    1636             :   !!              ier    - integer
    1637             :   !!                       ier = 0 normal and reliable termination of the
    1638             :   !!                               routine. it is assumed that the requested
    1639             :   !!                               accuracy has been achieved.
    1640             :   !!                       ier.gt.0 abnormal termination of the routine. the
    1641             :   !!                                estimates for result and error are less
    1642             :   !!                                reliable. it is assumed that the requested
    1643             :   !!                                accuracy has not been achieved.
    1644             :   !!
    1645             :   !!
    1646             :   !!***references  (none)
    1647             :   !!***routines called  sd1mach
    1648             :   !!***end prologue  dqk21
    1649             :   !!
    1650           0 :   subroutine dqk21(fx,fx_vars,a,b,result,abserr,resabs,resasc, ier)
    1651             : 
    1652             :     ! Arguments
    1653             :     interface
    1654             :       function fx(centr, vars)
    1655             :         use carma_precision_mod, only : f
    1656             :         use adgaquad_types_mod
    1657             :         real(kind=f), intent(in) :: centr
    1658             :         type(adgaquad_vars_type), intent(inout) :: vars
    1659             :         real(kind=f)             :: fx
    1660             :       end function fx
    1661             :     end interface
    1662             :     type(adgaquad_vars_type) :: fx_vars
    1663             :     real(kind=f) :: a
    1664             :     real(kind=f) :: b
    1665             :     real(kind=f) :: result
    1666             :     real(kind=f) :: abserr
    1667             :     real(kind=f) :: resabs
    1668             :     real(kind=f) :: resasc
    1669             :     integer :: ier
    1670             : 
    1671             :     ! Local declarations
    1672             :     real(kind=f) :: absc, centr, dabs, dhlgth, dmax1, dmin1
    1673             :     real(kind=f) :: epmach, fc, fsum, fval1, fval2, fv1(10), fv2(10), hlgth
    1674             :     real(kind=f) :: resg, resk, reskh, uflow, wg(5), wgk(11),xgk(11)
    1675             :     integer :: j,jtw,jtwm1
    1676             :  
    1677             :     !
    1678             :     !           the abscissae and weights are given for the interval (-1,1).
    1679             :     !           because of symmetry only the positive abscissae and their
    1680             :     !           corresponding weights are given.
    1681             :     !
    1682             :     !           xgk    - abscissae of the 21-point kronrod rule
    1683             :     !                    xgk(2), xgk(4), ...  abscissae of the 10-point
    1684             :     !                    gauss rule
    1685             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    1686             :     !                    added to the 10-point gauss rule
    1687             :     !
    1688             :     !           wgk    - weights of the 21-point kronrod rule
    1689             :     !
    1690             :     !           wg     - weights of the 10-point gauss rule
    1691             :     !
    1692             :     !
    1693             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    1694             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    1695             :     ! bell labs, nov. 1981.
    1696             :     !
    1697             :     data wg  (  1) / 0.066671344308688137593568809893332_f /
    1698             :     data wg  (  2) / 0.149451349150580593145776339657697_f /
    1699             :     data wg  (  3) / 0.219086362515982043995534934228163_f /
    1700             :     data wg  (  4) / 0.269266719309996355091226921569469_f /
    1701             :     data wg  (  5) / 0.295524224714752870173892994651338_f /
    1702             : 
    1703             :     data xgk (  1) / 0.995657163025808080735527280689003_f /
    1704             :     data xgk (  2) / 0.973906528517171720077964012084452_f /
    1705             :     data xgk (  3) / 0.930157491355708226001207180059508_f /
    1706             :     data xgk (  4) / 0.865063366688984510732096688423493_f /
    1707             :     data xgk (  5) / 0.780817726586416897063717578345042_f /
    1708             :     data xgk (  6) / 0.679409568299024406234327365114874_f /
    1709             :     data xgk (  7) / 0.562757134668604683339000099272694_f /
    1710             :     data xgk (  8) / 0.433395394129247190799265943165784_f /
    1711             :     data xgk (  9) / 0.294392862701460198131126603103866_f /
    1712             :     data xgk ( 10) / 0.148874338981631210884826001129720_f /
    1713             :     data xgk ( 11) / 0.000000000000000000000000000000000_f /
    1714             : 
    1715             :     data wgk (  1) / 0.011694638867371874278064396062192_f /
    1716             :     data wgk (  2) / 0.032558162307964727478818972459390_f /
    1717             :     data wgk (  3) / 0.054755896574351996031381300244580_f /
    1718             :     data wgk (  4) / 0.075039674810919952767043140916190_f /
    1719             :     data wgk (  5) / 0.093125454583697605535065465083366_f /
    1720             :     data wgk (  6) / 0.109387158802297641899210590325805_f /
    1721             :     data wgk (  7) / 0.123491976262065851077958109831074_f /
    1722             :     data wgk (  8) / 0.134709217311473325928054001771707_f /
    1723             :     data wgk (  9) / 0.142775938577060080797094273138717_f /
    1724             :     data wgk ( 10) / 0.147739104901338491374841515972068_f /
    1725             :     data wgk ( 11) / 0.149445554002916905664936468389821_f /
    1726             : 
    1727             :     !
    1728             :     !           list of major variables
    1729             :     !           -----------------------
    1730             :     !
    1731             :     !           centr  - mid point of the interval
    1732             :     !           hlgth  - half-length of the interval
    1733             :     !           absc   - abscissa
    1734             :     !           fval*  - function value
    1735             :     !           resg   - result of the 10-point gauss formula
    1736             :     !           resk   - result of the 21-point kronrod formula
    1737             :     !           reskh  - approximation to the mean value of f over (a,b),
    1738             :     !                    i.e. to i/(b-a)
    1739             :     !
    1740             :     !
    1741             :     !           machine dependent constants
    1742             :     !           ---------------------------
    1743             :     !
    1744             :     !           epmach is the largest relative spacing.
    1745             :     !           uflow is the smallest positive magnitude.
    1746             :     !
    1747             :     !***first executable statement  dqk21
    1748             :     !epmach = d1mach(4)
    1749             :     !uflow = d1mach(1)
    1750             : 
    1751           0 :     call sd1mach(4,epmach,ier)
    1752           0 :     if(ier.eq.9) return
    1753           0 :     call sd1mach(1,uflow,ier)
    1754           0 :     if(ier.eq.9) return
    1755             : 
    1756           0 :     centr = 0.5_f*(a+b)
    1757           0 :     hlgth = 0.5_f*(b-a)
    1758           0 :     dhlgth = dabs(hlgth)
    1759             :     !
    1760             :     !           compute the 21-point kronrod approximation to
    1761             :     !           the integral, and estimate the absolute error.
    1762             :     !
    1763           0 :     resg = 0.0_f
    1764           0 :     fc = fx(centr, fx_vars)
    1765           0 :     resk = wgk(11)*fc
    1766           0 :     resabs = dabs(resk)
    1767           0 :     do j=1,5
    1768           0 :       jtw = 2*j
    1769           0 :       absc = hlgth*xgk(jtw)
    1770           0 :       fval1 = fx(centr-absc, fx_vars)
    1771           0 :       fval2 = fx(centr+absc, fx_vars)
    1772           0 :       fv1(jtw) = fval1
    1773           0 :       fv2(jtw) = fval2
    1774           0 :       fsum = fval1+fval2
    1775           0 :       resg = resg+wg(j)*fsum
    1776           0 :       resk = resk+wgk(jtw)*fsum
    1777           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    1778             :     end do
    1779           0 :     do j = 1,5
    1780           0 :       jtwm1 = 2*j-1
    1781           0 :       absc = hlgth*xgk(jtwm1)
    1782           0 :       fval1 = fx(centr-absc, fx_vars)
    1783           0 :       fval2 = fx(centr+absc, fx_vars)
    1784           0 :       fv1(jtwm1) = fval1
    1785           0 :       fv2(jtwm1) = fval2
    1786           0 :       fsum = fval1+fval2
    1787           0 :       resk = resk+wgk(jtwm1)*fsum
    1788           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    1789             :     end do
    1790           0 :     reskh = resk*0.5_f
    1791           0 :     resasc = wgk(11)*dabs(fc-reskh)
    1792           0 :     do j=1,10
    1793           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    1794             :     end do
    1795           0 :     result = resk*hlgth
    1796           0 :     resabs = resabs*dhlgth
    1797           0 :     resasc = resasc*dhlgth
    1798           0 :     abserr = dabs((resk-resg)*hlgth)
    1799           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0.0_f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    1800           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    1801             : 999 return
    1802             :   end subroutine dqk21
    1803             : 
    1804             :   !!***begin prologue  dqk31
    1805             :   !!***date written   800101   (yymmdd)
    1806             :   !!***revision date  130519   (yymmdd)
    1807             :   !!***category no.  h2a1a2
    1808             :   !!***keywords  31-point gauss-kronrod rules
    1809             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    1810             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    1811             :   !!***purpose  to compute i = integral of f over (a,b) with error
    1812             :   !!                           estimate
    1813             :   !!                       j = integral of abs(f) over (a,b)
    1814             :   !!***description
    1815             :   !!
    1816             :   !!           integration rules
    1817             :   !!           standard fortran subroutine
    1818             :   !!           double precision version
    1819             :   !!
    1820             :   !!           parameters
    1821             :   !!            on entry
    1822             :   !!              fx     - double precision
    1823             :   !!                       function subprogram defining the integrand
    1824             :   !!                       function f(x). the actual name for f needs to be
    1825             :   !!                       declared e x t e r n a l in the calling program.
    1826             :   !!
    1827             :   !!            fx_vars- structure containing variables need for integration
    1828             :   !!                     specific to fractal meanfield scattering code!
    1829             :   !!
    1830             :   !!              a      - double precision
    1831             :   !!                       lower limit of integration
    1832             :   !!
    1833             :   !!              b      - double precision
    1834             :   !!                       upper limit of integration
    1835             :   !!
    1836             :   !!            on return
    1837             :   !!              result - double precision
    1838             :   !!                       approximation to the integral i
    1839             :   !!                       result is computed by applying the 31-point
    1840             :   !!                       gauss-kronrod rule (resk), obtained by optimal
    1841             :   !!                       addition of abscissae to the 15-point gauss
    1842             :   !!                       rule (resg).
    1843             :   !!
    1844             :   !!              abserr - double precison
    1845             :   !!                       estimate of the modulus of the modulus,
    1846             :   !!                       which should not exceed abs(i-result)
    1847             :   !!
    1848             :   !!              resabs - double precision
    1849             :   !!                       approximation to the integral j
    1850             :   !!
    1851             :   !!              resasc - double precision
    1852             :   !!                       approximation to the integral of abs(f-i/(b-a))
    1853             :   !!                       over (a,b)
    1854             :   !!
    1855             :   !!              ier    - integer
    1856             :   !!                       ier = 0 normal and reliable termination of the
    1857             :   !!                               routine. it is assumed that the requested
    1858             :   !!                               accuracy has been achieved.
    1859             :   !!                       ier.gt.0 abnormal termination of the routine. the
    1860             :   !!                                estimates for result and error are less
    1861             :   !!                                reliable. it is assumed that the requested
    1862             :   !!                                accuracy has not been achieved.
    1863             :   !!
    1864             :   !!
    1865             :   !!***references  (none)
    1866             :   !!***routines called  sd1mach
    1867             :   !!***end prologue  dqk31
    1868           0 :   subroutine dqk31(fx,fx_vars,a,b,result,abserr,resabs,resasc, ier)
    1869             : 
    1870             :     ! Arguments
    1871             :     interface
    1872             :       function fx(centr, vars)
    1873             :         use carma_precision_mod, only : f
    1874             :         use adgaquad_types_mod
    1875             :         real(kind=f), intent(in) :: centr
    1876             :         type(adgaquad_vars_type), intent(inout) :: vars
    1877             :         real(kind=f)             :: fx
    1878             :       end function fx
    1879             :     end interface
    1880             :     type(adgaquad_vars_type) :: fx_vars
    1881             :     real(kind=f) :: a
    1882             :     real(kind=f) :: b
    1883             :     real(kind=f) :: result
    1884             :     real(kind=f) :: abserr
    1885             :     real(kind=f) :: resabs
    1886             :     real(kind=f) :: resasc
    1887             :     integer :: ier
    1888             : 
    1889             :     ! Local declarations
    1890             :     real(kind=f) :: absc, centr, dabs, dhlgth, dmax1, dmin1
    1891             :     real(kind=f) :: epmach, fc, fsum, fval1, fval2, fv1(15), fv2(15), hlgth
    1892             :     real(kind=f) :: resg, resk, reskh, uflow, wg(8), wgk(16), xgk(16)
    1893             :     integer :: j,jtw,jtwm1
    1894             : 
    1895             :     !
    1896             :     !
    1897             :     !           the abscissae and weights are given for the interval (-1,1).
    1898             :     !           because of symmetry only the positive abscissae and their
    1899             :     !           corresponding weights are given.
    1900             :     !
    1901             :     !           xgk    - abscissae of the 31-point kronrod rule
    1902             :     !                    xgk(2), xgk(4), ...  abscissae of the 15-point
    1903             :     !                    gauss rule
    1904             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    1905             :     !                    added to the 15-point gauss rule
    1906             :     !
    1907             :     !           wgk    - weights of the 31-point kronrod rule
    1908             :     !
    1909             :     !           wg     - weights of the 15-point gauss rule
    1910             :     !
    1911             :     !
    1912             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    1913             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    1914             :     ! bell labs, nov. 1981.
    1915             :     !
    1916             :     data wg  (  1) / 0.030753241996117268354628393577204_f /
    1917             :     data wg  (  2) / 0.070366047488108124709267416450667_f /
    1918             :     data wg  (  3) / 0.107159220467171935011869546685869_f /
    1919             :     data wg  (  4) / 0.139570677926154314447804794511028_f /
    1920             :     data wg  (  5) / 0.166269205816993933553200860481209_f /
    1921             :     data wg  (  6) / 0.186161000015562211026800561866423_f /
    1922             :     data wg  (  7) / 0.198431485327111576456118326443839_f /
    1923             :     data wg  (  8) / 0.202578241925561272880620199967519_f /
    1924             : 
    1925             :     data xgk (  1) / 0.998002298693397060285172840152271_f /
    1926             :     data xgk (  2) / 0.987992518020485428489565718586613_f /
    1927             :     data xgk (  3) / 0.967739075679139134257347978784337_f /
    1928             :     data xgk (  4) / 0.937273392400705904307758947710209_f /
    1929             :     data xgk (  5) / 0.897264532344081900882509656454496_f /
    1930             :     data xgk (  6) / 0.848206583410427216200648320774217_f /
    1931             :     data xgk (  7) / 0.790418501442465932967649294817947_f /
    1932             :     data xgk (  8) / 0.724417731360170047416186054613938_f /
    1933             :     data xgk (  9) / 0.650996741297416970533735895313275_f /
    1934             :     data xgk ( 10) / 0.570972172608538847537226737253911_f /
    1935             :     data xgk ( 11) / 0.485081863640239680693655740232351_f /
    1936             :     data xgk ( 12) / 0.394151347077563369897207370981045_f /
    1937             :     data xgk ( 13) / 0.299180007153168812166780024266389_f /
    1938             :     data xgk ( 14) / 0.201194093997434522300628303394596_f /
    1939             :     data xgk ( 15) / 0.101142066918717499027074231447392_f /
    1940             :     data xgk ( 16) / 0.000000000000000000000000000000000_f /
    1941             : 
    1942             :     data wgk (  1) / 0.005377479872923348987792051430128_f /
    1943             :     data wgk (  2) / 0.015007947329316122538374763075807_f /
    1944             :     data wgk (  3) / 0.025460847326715320186874001019653_f /
    1945             :     data wgk (  4) / 0.035346360791375846222037948478360_f /
    1946             :     data wgk (  5) / 0.044589751324764876608227299373280_f /
    1947             :     data wgk (  6) / 0.053481524690928087265343147239430_f /
    1948             :     data wgk (  7) / 0.062009567800670640285139230960803_f /
    1949             :     data wgk (  8) / 0.069854121318728258709520077099147_f /
    1950             :     data wgk (  9) / 0.076849680757720378894432777482659_f /
    1951             :     data wgk ( 10) / 0.083080502823133021038289247286104_f /
    1952             :     data wgk ( 11) / 0.088564443056211770647275443693774_f /
    1953             :     data wgk ( 12) / 0.093126598170825321225486872747346_f /
    1954             :     data wgk ( 13) / 0.096642726983623678505179907627589_f /
    1955             :     data wgk ( 14) / 0.099173598721791959332393173484603_f /
    1956             :     data wgk ( 15) / 0.100769845523875595044946662617570_f /
    1957             :     data wgk ( 16) / 0.101330007014791549017374792767493_f /
    1958             :     !
    1959             :     !
    1960             :     !           list of major variables
    1961             :     !           -----------------------
    1962             :     !           centr  - mid point of the interval
    1963             :     !           hlgth  - half-length of the interval
    1964             :     !           absc   - abscissa
    1965             :     !           fval*  - function value
    1966             :     !           resg   - result of the 15-point gauss formula
    1967             :     !           resk   - result of the 31-point kronrod formula
    1968             :     !           reskh  - approximation to the mean value of f over (a,b),
    1969             :     !                    i.e. to i/(b-a)
    1970             :     !
    1971             :     !           machine dependent constants
    1972             :     !           ---------------------------
    1973             :     !           epmach is the largest relative spacing.
    1974             :     !           uflow is the smallest positive magnitude.
    1975             :     !***first executable statement  dqk31
    1976           0 :     call sd1mach(4,epmach,ier)
    1977           0 :     if(ier.eq.9) return
    1978           0 :     call sd1mach(1,uflow,ier)
    1979           0 :     if(ier.eq.9) return
    1980             : 
    1981             :     !epmach = d1mach(4)
    1982             :     !uflow = d1mach(1)
    1983             : 
    1984           0 :     centr = 0.5_f*(a+b)
    1985           0 :     hlgth = 0.5_f*(b-a)
    1986           0 :     dhlgth = dabs(hlgth)
    1987             :     !
    1988             :     !           compute the 31-point kronrod approximation to
    1989             :     !           the integral, and estimate the absolute error.
    1990             :     !
    1991           0 :     fc = fx(centr, fx_vars)
    1992           0 :     resg = wg(8)*fc
    1993           0 :     resk = wgk(16)*fc
    1994           0 :     resabs = dabs(resk)
    1995           0 :     do j=1,7
    1996           0 :       jtw = j*2
    1997           0 :       absc = hlgth*xgk(jtw)
    1998           0 :       fval1 = fx(centr-absc, fx_vars)
    1999           0 :       fval2 = fx(centr+absc, fx_vars)
    2000           0 :       fv1(jtw) = fval1
    2001           0 :       fv2(jtw) = fval2
    2002           0 :       fsum = fval1+fval2
    2003           0 :       resg = resg+wg(j)*fsum
    2004           0 :       resk = resk+wgk(jtw)*fsum
    2005           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    2006             :     end do
    2007           0 :     do j = 1,8
    2008           0 :       jtwm1 = j*2-1
    2009           0 :       absc = hlgth*xgk(jtwm1)
    2010           0 :       fval1 = fx(centr-absc, fx_vars)
    2011           0 :       fval2 = fx(centr+absc, fx_vars)
    2012           0 :       fv1(jtwm1) = fval1
    2013           0 :       fv2(jtwm1) = fval2
    2014           0 :       fsum = fval1+fval2
    2015           0 :       resk = resk+wgk(jtwm1)*fsum
    2016           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    2017             :     end do
    2018           0 :     reskh = resk*0.5_f
    2019           0 :     resasc = wgk(16)*dabs(fc-reskh)
    2020           0 :     do j=1,15
    2021           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    2022             :     end do
    2023           0 :     result = resk*hlgth
    2024           0 :     resabs = resabs*dhlgth
    2025           0 :     resasc = resasc*dhlgth
    2026           0 :     abserr = dabs((resk-resg)*hlgth)
    2027           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0.0_f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    2028           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    2029             : 999 return
    2030             :   end subroutine dqk31
    2031             : 
    2032             : 
    2033             : 
    2034             :   !!***begin prologue  dqk41
    2035             :   !!***date written   800101   (yymmdd)
    2036             :   !!***revision date  130319   (yymmdd)
    2037             :   !!***category no.  h2a1a2
    2038             :   !!***keywords  41-point gauss-kronrod rules
    2039             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    2040             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    2041             :   !!***purpose  to compute i = integral of f over (a,b), with error
    2042             :   !!                           estimate
    2043             :   !!                       j = integral of abs(f) over (a,b)
    2044             :   !!***description
    2045             :   !!
    2046             :   !!           integration rules
    2047             :   !!           standard fortran subroutine
    2048             :   !!           double precision version
    2049             :   !!
    2050             :   !!           parameters
    2051             :   !!            on entry
    2052             :   !!              fx     - double precision
    2053             :   !!                       function subprogram defining the integrand
    2054             :   !!                       function f(x). the actual name for f needs to be
    2055             :   !!                       declared e x t e r n a l in the calling program.
    2056             :   !!
    2057             :   !!            fx_vars- structure containing variables need for integration
    2058             :   !!                     specific to fractal meanfield scattering code
    2059             :   !!
    2060             :   !!              a      - double precision
    2061             :   !!                       lower limit of integration
    2062             :   !!
    2063             :   !!              b      - double precision
    2064             :   !!                       upper limit of integration
    2065             :   !!
    2066             :   !!            on return
    2067             :   !!              result - double precision
    2068             :   !!                       approximation to the integral i
    2069             :   !!                       result is computed by applying the 41-point
    2070             :   !!                       gauss-kronrod rule (resk) obtained by optimal
    2071             :   !!                       addition of abscissae to the 20-point gauss
    2072             :   !!                       rule (resg).
    2073             :   !!
    2074             :   !!              abserr - double precision
    2075             :   !!                       estimate of the modulus of the absolute error,
    2076             :   !!                       which should not exceed abs(i-result)
    2077             :   !!
    2078             :   !!              resabs - double precision
    2079             :   !!                       approximation to the integral j
    2080             :   !!
    2081             :   !!              resasc - double precision
    2082             :   !!                       approximation to the integal of abs(f-i/(b-a))
    2083             :   !!                       over (a,b)
    2084             :   !!
    2085             :   !!              ier    - integer
    2086             :   !!                       ier = 0 normal and reliable termination of the
    2087             :   !!                               routine. it is assumed that the requested
    2088             :   !!                               accuracy has been achieved.
    2089             :   !!                       ier.gt.0 abnormal termination of the routine. the
    2090             :   !!                                estimates for result and error are less
    2091             :   !!                                reliable. it is assumed that the requested
    2092             :   !!                                accuracy has not been achieved.
    2093             :   !!
    2094             :   !!
    2095             :   !!***references  (none)
    2096             :   !!***routines called  sd1mach
    2097             :   !!***end prologue  dqk41
    2098             :   !!
    2099           0 :   subroutine dqk41(fx,fx_vars,a,b,result,abserr,resabs,resasc, ier)
    2100             : 
    2101             :     ! Arguments
    2102             :     interface
    2103             :       function fx(centr, vars)
    2104             :         use carma_precision_mod, only : f
    2105             :         use adgaquad_types_mod
    2106             :         real(kind=f), intent(in) :: centr
    2107             :         type(adgaquad_vars_type), intent(inout) :: vars
    2108             :         real(kind=f)             :: fx
    2109             :       end function fx
    2110             :     end interface
    2111             :     type(adgaquad_vars_type) :: fx_vars
    2112             :     real(kind=f) :: a
    2113             :     real(kind=f) :: b
    2114             :     real(kind=f) :: result
    2115             :     real(kind=f) :: abserr
    2116             :     real(kind=f) :: resabs
    2117             :     real(kind=f) :: resasc
    2118             :     integer :: ier
    2119             : 
    2120             :     ! Local declarations
    2121             :     real(kind=f) :: absc, centr, dabs, dhlgth, dmax1, dmin1
    2122             :     real(kind=f) :: epmach, fc, fsum, fval1, fval2, fv1(20), fv2(20), hlgth
    2123             :     real(kind=f) :: resg, resk, reskh, uflow, wg(10), wgk(21), xgk(21)
    2124             :     integer :: j, jtw, jtwm1
    2125             : 
    2126             :     !
    2127             :     !           the abscissae and weights are given for the interval (-1,1).
    2128             :     !           because of symmetry only the positive abscissae and their
    2129             :     !           corresponding weights are given.
    2130             :     !
    2131             :     !           xgk    - abscissae of the 41-point gauss-kronrod rule
    2132             :     !                    xgk(2), xgk(4), ...  abscissae of the 20-point
    2133             :     !                    gauss rule
    2134             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    2135             :     !                    added to the 20-point gauss rule
    2136             :     !
    2137             :     !           wgk    - weights of the 41-point gauss-kronrod rule
    2138             :     !
    2139             :     !           wg     - weights of the 20-point gauss rule
    2140             :     !
    2141             :     !
    2142             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    2143             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    2144             :     ! bell labs, nov. 1981.
    2145             :     !
    2146             :     data wg  (  1) / 0.017614007139152118311861962351853_f /
    2147             :     data wg  (  2) / 0.040601429800386941331039952274932_f /
    2148             :     data wg  (  3) / 0.062672048334109063569506535187042_f /
    2149             :     data wg  (  4) / 0.083276741576704748724758143222046_f /
    2150             :     data wg  (  5) / 0.101930119817240435036750135480350_f /
    2151             :     data wg  (  6) / 0.118194531961518417312377377711382_f /
    2152             :     data wg  (  7) / 0.131688638449176626898494499748163_f /
    2153             :     data wg  (  8) / 0.142096109318382051329298325067165_f /
    2154             :     data wg  (  9) / 0.149172986472603746787828737001969_f /
    2155             :     data wg  ( 10) / 0.152753387130725850698084331955098_f /
    2156             : 
    2157             :     data xgk (  1) / 0.998859031588277663838315576545863_f /
    2158             :     data xgk (  2) / 0.993128599185094924786122388471320_f /
    2159             :     data xgk (  3) / 0.981507877450250259193342994720217_f /
    2160             :     data xgk (  4) / 0.963971927277913791267666131197277_f /
    2161             :     data xgk (  5) / 0.940822633831754753519982722212443_f /
    2162             :     data xgk (  6) / 0.912234428251325905867752441203298_f /
    2163             :     data xgk (  7) / 0.878276811252281976077442995113078_f /
    2164             :     data xgk (  8) / 0.839116971822218823394529061701521_f /
    2165             :     data xgk (  9) / 0.795041428837551198350638833272788_f /
    2166             :     data xgk ( 10) / 0.746331906460150792614305070355642_f /
    2167             :     data xgk ( 11) / 0.693237656334751384805490711845932_f /
    2168             :     data xgk ( 12) / 0.636053680726515025452836696226286_f /
    2169             :     data xgk ( 13) / 0.575140446819710315342946036586425_f /
    2170             :     data xgk ( 14) / 0.510867001950827098004364050955251_f/
    2171             :     data xgk ( 15) / 0.443593175238725103199992213492640_f /
    2172             :     data xgk ( 16) / 0.373706088715419560672548177024927_f /
    2173             :     data xgk ( 17) / 0.301627868114913004320555356858592_f /
    2174             :     data xgk ( 18) / 0.227785851141645078080496195368575_f /
    2175             :     data xgk ( 19) / 0.152605465240922675505220241022678_f /
    2176             :     data xgk ( 20) / 0.076526521133497333754640409398838_f /
    2177             :     data xgk ( 21) / 0.000000000000000000000000000000000_f /
    2178             : 
    2179             :     data wgk (  1) / 0.003073583718520531501218293246031_f /
    2180             :     data wgk (  2) / 0.008600269855642942198661787950102_f /
    2181             :     data wgk (  3) / 0.014626169256971252983787960308868_f /
    2182             :     data wgk (  4) / 0.020388373461266523598010231432755_f /
    2183             :     data wgk (  5) / 0.025882133604951158834505067096153_f /
    2184             :     data wgk (  6) / 0.031287306777032798958543119323801_f /
    2185             :     data wgk (  7) / 0.036600169758200798030557240707211_f /
    2186             :     data wgk (  8) / 0.041668873327973686263788305936895_f /
    2187             :     data wgk (  9) / 0.046434821867497674720231880926108_f /
    2188             :     data wgk ( 10) / 0.050944573923728691932707670050345_f /
    2189             :     data wgk ( 11) / 0.055195105348285994744832372419777_f /
    2190             :     data wgk ( 12) / 0.059111400880639572374967220648594_f /
    2191             :     data wgk ( 13) / 0.062653237554781168025870122174255_f /
    2192             :     data wgk ( 14) / 0.065834597133618422111563556969398_f /
    2193             :     data wgk ( 15) / 0.068648672928521619345623411885368_f /
    2194             :     data wgk ( 16) / 0.071054423553444068305790361723210_f /
    2195             :     data wgk ( 17) / 0.073030690332786667495189417658913_f /
    2196             :     data wgk ( 18) / 0.074582875400499188986581418362488_f /
    2197             :     data wgk ( 19) / 0.075704497684556674659542775376617_f /
    2198             :     data wgk ( 20) / 0.076377867672080736705502835038061_f /
    2199             :     data wgk ( 21) / 0.076600711917999656445049901530102_f /
    2200             :    
    2201             :     !
    2202             :     !           list of major variables
    2203             :     !           -----------------------
    2204             :     !
    2205             :     !           centr  - mid point of the interval
    2206             :     !           hlgth  - half-length of the interval
    2207             :     !           absc   - abscissa
    2208             :     !           fval*  - function value
    2209             :     !           resg   - result of the 20-point gauss formula
    2210             :     !           resk   - result of the 41-point kronrod formula
    2211             :     !           reskh  - approximation to mean value of f over (a,b), i.e.
    2212             :     !                    to i/(b-a)
    2213             :     !
    2214             :     !           machine dependent constants
    2215             :     !           ---------------------------
    2216             :     !
    2217             :     !           epmach is the largest relative spacing.
    2218             :     !           uflow is the smallest positive magnitude.
    2219             :     !
    2220             :     !***first executable statement  dqk41
    2221           0 :     call sd1mach(4,epmach,ier)
    2222           0 :     if(ier.eq.9) return
    2223           0 :     call sd1mach(1,uflow,ier)
    2224           0 :     if(ier.eq.9) return
    2225             : 
    2226             :     !epmach = d1mach(4)
    2227             :     !uflow = d1mach(1)
    2228             : 
    2229           0 :     centr = 0.5_f*(a+b)
    2230           0 :     hlgth = 0.5_f*(b-a)
    2231           0 :     dhlgth = dabs(hlgth)
    2232             :     !
    2233             :     !           compute the 41-point gauss-kronrod approximation to
    2234             :     !           the integral, and estimate the absolute error.
    2235             :     !
    2236           0 :     resg = 0.0_f
    2237           0 :     fc = fx(centr, fx_vars)
    2238           0 :     resk = wgk(21)*fc
    2239           0 :     resabs = dabs(resk)
    2240           0 :     do j=1,10
    2241           0 :       jtw = j*2
    2242           0 :       absc = hlgth*xgk(jtw)
    2243           0 :       fval1 = fx(centr-absc, fx_vars)
    2244           0 :       fval2 = fx(centr+absc, fx_vars)
    2245           0 :       fv1(jtw) = fval1
    2246           0 :       fv2(jtw) = fval2
    2247           0 :       fsum = fval1+fval2
    2248           0 :       resg = resg+wg(j)*fsum
    2249           0 :       resk = resk+wgk(jtw)*fsum
    2250           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    2251             :     end do
    2252           0 :     do j = 1,10
    2253           0 :       jtwm1 = j*2-1
    2254           0 :       absc = hlgth*xgk(jtwm1)
    2255           0 :       fval1 = fx(centr-absc, fx_vars)
    2256           0 :       fval2 = fx(centr+absc, fx_vars)
    2257           0 :       fv1(jtwm1) = fval1
    2258           0 :       fv2(jtwm1) = fval2
    2259           0 :       fsum = fval1+fval2
    2260           0 :       resk = resk+wgk(jtwm1)*fsum
    2261           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    2262             :     end do
    2263           0 :     reskh = resk*0.5_f
    2264           0 :     resasc = wgk(21)*dabs(fc-reskh)
    2265           0 :     do j=1,20
    2266           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    2267             :     end do
    2268           0 :     result = resk*hlgth
    2269           0 :     resabs = resabs*dhlgth
    2270           0 :     resasc = resasc*dhlgth
    2271           0 :     abserr = dabs((resk-resg)*hlgth)
    2272           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0._f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    2273           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    2274             : 999 return
    2275             :   end subroutine dqk41
    2276             : 
    2277             : 
    2278             :   !!***begin prologue  dqk51
    2279             :   !!***date written   800101   (yymmdd)
    2280             :   !!***revision date  130319   (yymmdd)
    2281             :   !!***category no.  h2a1a2
    2282             :   !!***keywords  51-point gauss-kronrod rules
    2283             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    2284             :   !!           de doncker,elise,appl. math & progr. div. - k.u.leuven
    2285             :   !!***purpose  to compute i = integral of f over (a,b) with error
    2286             :   !!                           estimate
    2287             :   !!                       j = integral of abs(f) over (a,b)
    2288             :   !!***description
    2289             :   !!
    2290             :   !!           integration rules
    2291             :   !!           standard fortran subroutine
    2292             :   !!           double precision version
    2293             :   !!
    2294             :   !!           parameters
    2295             :   !!            on entry
    2296             :   !!              fx     - double precision
    2297             :   !!                       function subroutine defining the integrand
    2298             :   !!                       function f(x). the actual name for f needs to be
    2299             :   !!                       declared e x t e r n a l in the calling program.
    2300             :   !!
    2301             :   !!              fx_vars- structure containing variables need for integration
    2302             :   !!                       specific to fractal meanfield scattering code
    2303             :   !!
    2304             :   !!              a      - double precision
    2305             :   !!                       lower limit of integration
    2306             :   !!
    2307             :   !!              b      - double precision
    2308             :   !!                       upper limit of integration
    2309             :   !!
    2310             :   !!            on return
    2311             :   !!              result - double precision
    2312             :   !!                       approximation to the integral i
    2313             :   !!                       result is computed by applying the 51-point
    2314             :   !!                       kronrod rule (resk) obtained by optimal addition
    2315             :   !!                       of abscissae to the 25-point gauss rule (resg).
    2316             :   !!
    2317             :   !!              abserr - double precision
    2318             :   !!                       estimate of the modulus of the absolute error,
    2319             :   !!                       which should not exceed abs(i-result)
    2320             :   !!
    2321             :   !!              resabs - double precision
    2322             :   !!                       approximation to the integral j
    2323             :   !!
    2324             :   !!              resasc - double precision
    2325             :   !!                       approximation to the integral of abs(f-i/(b-a))
    2326             :   !!                       over (a,b)
    2327             :   !!
    2328             :   !!              ier    - integer
    2329             :   !!                       ier = 0 normal and reliable termination of the
    2330             :   !!                               routine. it is assumed that the requested
    2331             :   !!                               accuracy has been achieved.
    2332             :   !!                       ier.gt.0 abnormal termination of the routine. the
    2333             :   !!                                estimates for result and error are less
    2334             :   !!                                reliable. it is assumed that the requested
    2335             :   !!                                accuracy has not been achieved.
    2336             :   !!
    2337             :   !!***references  (none)
    2338             :   !!***routines called  sd1mach
    2339             :   !!***end prologue  dqk51
    2340             :   !!
    2341           0 :   subroutine dqk51(fx,fx_vars,a,b,result,abserr,resabs,resasc,ier)
    2342             : 
    2343             :     ! Arguments
    2344             :     interface
    2345             :       function fx(centr, vars)
    2346             :         use carma_precision_mod, only : f
    2347             :         use adgaquad_types_mod
    2348             :         real(kind=f), intent(in) :: centr
    2349             :         type(adgaquad_vars_type), intent(inout) :: vars
    2350             :         real(kind=f)             :: fx
    2351             :       end function fx
    2352             :     end interface
    2353             :     type(adgaquad_vars_type) :: fx_vars
    2354             :     real(kind=f) :: a
    2355             :     real(kind=f) :: b
    2356             :     real(kind=f) :: result
    2357             :     real(kind=f) :: abserr
    2358             :     real(kind=f) :: resabs
    2359             :     real(kind=f) :: resasc
    2360             :     integer :: ier
    2361             : 
    2362             :     ! Local declarations
    2363             :     real(kind=f) ::  absc, centr, dabs, dhlgth, dmax1, dmin1
    2364             :     real(kind=f) ::  epmach, fc, fsum, fval1, fval2, fv1(25), fv2(25), hlgth
    2365             :     real(kind=f) ::  resg, resk, reskh, uflow, wg(13),wgk(26), xgk(26)
    2366             :     integer j,jtw,jtwm1
    2367             : 
    2368             :     !
    2369             :     !           the abscissae and weights are given for the interval (-1,1).
    2370             :     !           because of symmetry only the positive abscissae and their
    2371             :     !           corresponding weights are given.
    2372             :     !
    2373             :     !           xgk    - abscissae of the 51-point kronrod rule
    2374             :     !                    xgk(2), xgk(4), ...  abscissae of the 25-point
    2375             :     !                    gauss rule
    2376             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    2377             :     !                    added to the 25-point gauss rule
    2378             :     !
    2379             :     !           wgk    - weights of the 51-point kronrod rule
    2380             :     !
    2381             :     !           wg     - weights of the 25-point gauss rule
    2382             :     !
    2383             :     !
    2384             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    2385             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    2386             :     ! bell labs, nov. 1981.
    2387             :     !
    2388             :     data wg  (  1) / 0.011393798501026287947902964113235_f /
    2389             :     data wg  (  2) / 0.026354986615032137261901815295299_f /
    2390             :     data wg  (  3) / 0.040939156701306312655623487711646_f /
    2391             :     data wg  (  4) / 0.054904695975835191925936891540473_f /
    2392             :     data wg  (  5) / 0.068038333812356917207187185656708_f /
    2393             :     data wg  (  6) / 0.080140700335001018013234959669111_f /
    2394             :     data wg  (  7) / 0.091028261982963649811497220702892_f /
    2395             :     data wg  (  8) / 0.100535949067050644202206890392686_f /
    2396             :     data wg  (  9) / 0.108519624474263653116093957050117_f /
    2397             :     data wg  ( 10) / 0.114858259145711648339325545869556_f /
    2398             :     data wg  ( 11) / 0.119455763535784772228178126512901_f /
    2399             :     data wg  ( 12) / 0.122242442990310041688959518945852_f /
    2400             :     data wg  ( 13) / 0.123176053726715451203902873079050_f /
    2401             : 
    2402             :     data xgk (  1) / 0.999262104992609834193457486540341_f /
    2403             :     data xgk (  2) / 0.995556969790498097908784946893902_f /
    2404             :     data xgk (  3) / 0.988035794534077247637331014577406_f /
    2405             :     data xgk (  4) / 0.976663921459517511498315386479594_f /
    2406             :     data xgk (  5) / 0.961614986425842512418130033660167_f /
    2407             :     data xgk (  6) / 0.942974571228974339414011169658471_f /
    2408             :     data xgk (  7) / 0.920747115281701561746346084546331_f /
    2409             :     data xgk (  8) / 0.894991997878275368851042006782805_f /
    2410             :     data xgk (  9) / 0.865847065293275595448996969588340_f /
    2411             :     data xgk ( 10) / 0.833442628760834001421021108693570_f /
    2412             :     data xgk ( 11) / 0.797873797998500059410410904994307_f /
    2413             :     data xgk ( 12) / 0.759259263037357630577282865204361_f /
    2414             :     data xgk ( 13) / 0.717766406813084388186654079773298_f /
    2415             :     data xgk ( 14) / 0.673566368473468364485120633247622_f /
    2416             :     data xgk ( 15) / 0.626810099010317412788122681624518_f /
    2417             :     data xgk ( 16) / 0.577662930241222967723689841612654_f /
    2418             :     data xgk ( 17) / 0.526325284334719182599623778158010_f /
    2419             :     data xgk ( 18) / 0.473002731445714960522182115009192_f /
    2420             :     data xgk ( 19) / 0.417885382193037748851814394594572_f /
    2421             :     data xgk ( 20) / 0.361172305809387837735821730127641_f /
    2422             :     data xgk ( 21) / 0.303089538931107830167478909980339_f /
    2423             :     data xgk ( 22) / 0.243866883720988432045190362797452_f /
    2424             :     data xgk ( 23) / 0.183718939421048892015969888759528_f /
    2425             :     data xgk ( 24) / 0.122864692610710396387359818808037_f /
    2426             :     data xgk ( 25) / 0.061544483005685078886546392366797_f /
    2427             :     data xgk ( 26) / 0.000000000000000000000000000000000_f /
    2428             : 
    2429             :     data wgk (  1) / 0.001987383892330315926507851882843_f /
    2430             :     data wgk (  2) / 0.005561932135356713758040236901066_f /
    2431             :     data wgk (  3) / 0.009473973386174151607207710523655_f /
    2432             :     data wgk (  4) / 0.013236229195571674813656405846976_f /
    2433             :     data wgk (  5) / 0.016847817709128298231516667536336_f /
    2434             :     data wgk (  6) / 0.020435371145882835456568292235939_f /
    2435             :     data wgk (  7) / 0.024009945606953216220092489164881_f /
    2436             :     data wgk (  8) / 0.027475317587851737802948455517811_f /
    2437             :     data wgk (  9) / 0.030792300167387488891109020215229_f /
    2438             :     data wgk ( 10) / 0.034002130274329337836748795229551_f /
    2439             :     data wgk ( 11) / 0.037116271483415543560330625367620_f /
    2440             :     data wgk ( 12) / 0.040083825504032382074839284467076_f /
    2441             :     data wgk ( 13) / 0.042872845020170049476895792439495_f /
    2442             :     data wgk ( 14) / 0.045502913049921788909870584752660_f /
    2443             :     data wgk ( 15) / 0.047982537138836713906392255756915_f /
    2444             :     data wgk ( 16) / 0.050277679080715671963325259433440_f /
    2445             :     data wgk ( 17) / 0.052362885806407475864366712137873_f /
    2446             :     data wgk ( 18) / 0.054251129888545490144543370459876_f /
    2447             :     data wgk ( 19) / 0.055950811220412317308240686382747_f /
    2448             :     data wgk ( 20) / 0.057437116361567832853582693939506_f /
    2449             :     data wgk ( 21) / 0.058689680022394207961974175856788_f /
    2450             :     data wgk ( 22) / 0.059720340324174059979099291932562_f /
    2451             :     data wgk ( 23) / 0.060539455376045862945360267517565_f /
    2452             :     data wgk ( 24) / 0.061128509717053048305859030416293_f /
    2453             :     data wgk ( 25) / 0.061471189871425316661544131965264_f /
    2454             :     !       note: wgk (26) was calculated from the values of wgk(1..25)
    2455             :     data wgk ( 26) / 0.061580818067832935078759824240066_f /
    2456             :    
    2457             :     !
    2458             :     !           list of major variables
    2459             :     !           -----------------------
    2460             :     !
    2461             :     !           centr  - mid point of the interval
    2462             :     !           hlgth  - half-length of the interval
    2463             :     !           absc   - abscissa
    2464             :     !           fval*  - function value
    2465             :     !           resg   - result of the 25-point gauss formula
    2466             :     !           resk   - result of the 51-point kronrod formula
    2467             :     !           reskh  - approximation to the mean value of f over (a,b),
    2468             :     !                    i.e. to i/(b-a)
    2469             :     !
    2470             :     !           machine dependent constants
    2471             :     !           ---------------------------
    2472             :     !
    2473             :     !           epmach is the largest relative spacing.
    2474             :     !           uflow is the smallest positive magnitude.
    2475             :     !
    2476             :     !***first executable statement  dqk51
    2477           0 :     call sd1mach(4,epmach,ier)
    2478           0 :     if(ier.eq.9) return
    2479           0 :     call sd1mach(1,uflow,ier)
    2480           0 :     if(ier.eq.9) return
    2481             : 
    2482             :     !epmach = d1mach(4)
    2483             :     !uflow = d1mach(1)
    2484             : 
    2485           0 :     centr = 0.5_f*(a+b)
    2486           0 :     hlgth = 0.5_f*(b-a)
    2487           0 :     dhlgth = dabs(hlgth)
    2488             :     !
    2489             :     !           compute the 51-point kronrod approximation to
    2490             :     !           the integral, and estimate the absolute error.
    2491             :     !
    2492           0 :     fc = fx(centr, fx_vars)
    2493           0 :     resg = wg(13)*fc
    2494           0 :     resk = wgk(26)*fc
    2495           0 :     resabs = dabs(resk)
    2496           0 :     do j=1,12
    2497           0 :       jtw = j*2
    2498           0 :       absc = hlgth*xgk(jtw)
    2499           0 :       fval1 = fx(centr-absc, fx_vars)
    2500           0 :       fval2 = fx(centr+absc, fx_vars)
    2501           0 :       fv1(jtw) = fval1
    2502           0 :       fv2(jtw) = fval2
    2503           0 :       fsum = fval1+fval2
    2504           0 :       resg = resg+wg(j)*fsum
    2505           0 :       resk = resk+wgk(jtw)*fsum
    2506           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    2507             :     end do
    2508           0 :     do j = 1,13
    2509           0 :       jtwm1 = j*2-1
    2510           0 :       absc = hlgth*xgk(jtwm1)
    2511           0 :       fval1 = fx(centr-absc, fx_vars)
    2512           0 :       fval2 = fx(centr+absc, fx_vars)
    2513           0 :       fv1(jtwm1) = fval1
    2514           0 :       fv2(jtwm1) = fval2
    2515           0 :       fsum = fval1+fval2
    2516           0 :       resk = resk+wgk(jtwm1)*fsum
    2517           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    2518             :     end do
    2519           0 :     reskh = resk*0.5_f
    2520           0 :     resasc = wgk(26)*dabs(fc-reskh)
    2521           0 :     do j=1,25
    2522           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    2523             :     end do
    2524           0 :     result = resk*hlgth
    2525           0 :     resabs = resabs*dhlgth
    2526           0 :     resasc = resasc*dhlgth
    2527           0 :     abserr = dabs((resk-resg)*hlgth)
    2528           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0.0_f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    2529           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    2530             : 999 return
    2531             :   end subroutine dqk51
    2532             : 
    2533             :   !!***begin prologue  dqk61
    2534             :   !!***date written   800101   (yymmdd)
    2535             :   !!***revision date  130319   (yymmdd)
    2536             :   !!***category no.  h2a1a2
    2537             :   !!***keywords  61-point gauss-kronrod rules
    2538             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    2539             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    2540             :   !!***purpose  to compute i = integral of f over (a,b) with error
    2541             :   !!                           estimate
    2542             :   !!                       j = integral of dabs(f) over (a,b)
    2543             :   !!***description
    2544             :   !!
    2545             :   !!        integration rule
    2546             :   !!        standard fortran subroutine
    2547             :   !!        double precision version
    2548             :   !!
    2549             :   !!
    2550             :   !!        parameters
    2551             :   !!         on entry
    2552             :   !!           fx     - double precision
    2553             :   !!                    function subprogram defining the integrand
    2554             :   !!                    function f(x). the actual name for f needs to be
    2555             :   !!                    declared e x t e r n a l in the calling program.
    2556             :   !!
    2557             :   !!            fx_vars- structure containing variables need for integration
    2558             :   !!                     specific to fractal meanfield scattering code
    2559             :   !!
    2560             :   !!           a      - double precision
    2561             :   !!                    lower limit of integration
    2562             :   !!
    2563             :   !!           b      - double precision
    2564             :   !!                    upper limit of integration
    2565             :   !!
    2566             :   !!         on return
    2567             :   !!           result - double precision
    2568             :   !!                    approximation to the integral i
    2569             :   !!                    result is computed by applying the 61-point
    2570             :   !!                    kronrod rule (resk) obtained by optimal addition of
    2571             :   !!                    abscissae to the 30-point gauss rule (resg).
    2572             :   !!
    2573             :   !!           abserr - double precision
    2574             :   !!                    estimate of the modulus of the absolute error,
    2575             :   !!                    which should equal or exceed dabs(i-result)
    2576             :   !!
    2577             :   !!           resabs - double precision
    2578             :   !!                    approximation to the integral j
    2579             :   !!
    2580             :   !!           resasc - double precision
    2581             :   !!                    approximation to the integral of dabs(f-i/(b-a))
    2582             :   !!
    2583             :   !!              ier    - integer
    2584             :   !!                       ier = 0 normal and reliable termination of the
    2585             :   !!                               routine. it is assumed that the requested
    2586             :   !!                               accuracy has been achieved.
    2587             :   !!                       ier.gt.0 abnormal termination of the routine. the
    2588             :   !!                                estimates for result and error are less
    2589             :   !!                                reliable. it is assumed that the requested
    2590             :   !!                                accuracy has not been achieved.
    2591             :   !!
    2592             :   !!
    2593             :   !!***references  (none)
    2594             :   !!***routines called  sd1mach
    2595             :   !!***end prologue  dqk61
    2596             :   !!
    2597           0 :   subroutine dqk61(fx,fx_vars,a,b,result,abserr,resabs,resasc, ier)
    2598             : 
    2599             :     ! Arguments
    2600             :     interface
    2601             :       function fx(centr, vars)
    2602             :         use carma_precision_mod, only : f
    2603             :         use adgaquad_types_mod
    2604             :         real(kind=f), intent(in) :: centr
    2605             :         type(adgaquad_vars_type), intent(inout) :: vars
    2606             :         real(kind=f)             :: fx
    2607             :       end function fx
    2608             :     end interface
    2609             :     type(adgaquad_vars_type) :: fx_vars
    2610             :     real(kind=f) :: a
    2611             :     real(kind=f) :: b
    2612             :     real(kind=f) :: result
    2613             :     real(kind=f) :: abserr
    2614             :     real(kind=f) :: resabs
    2615             :     real(kind=f) :: resasc
    2616             :     integer :: ier
    2617             : 
    2618             :     ! Local declartions
    2619             :     real(kind=f) :: dabsc, centr, dabs, dhlgth, dmax1, dmin1
    2620             :     real(kind=f) :: epmach, fc, fsum, fval1, fval2, fv1(30), fv2(30), hlgth
    2621             :     real(kind=f) :: resg, resk, reskh, uflow, wg(15) ,wgk(31), xgk(31)
    2622             :     integer :: j, jtw, jtwm1
    2623             : 
    2624             :     !
    2625             :     !           the abscissae and weights are given for the
    2626             :     !           interval (-1,1). because of symmetry only the positive
    2627             :     !           abscissae and their corresponding weights are given.
    2628             :     !
    2629             :     !           xgk   - abscissae of the 61-point kronrod rule
    2630             :     !                   xgk(2), xgk(4)  ... abscissae of the 30-point
    2631             :     !                   gauss rule
    2632             :     !                   xgk(1), xgk(3)  ... optimally added abscissae
    2633             :     !                   to the 30-point gauss rule
    2634             :     !
    2635             :     !           wgk   - weights of the 61-point kronrod rule
    2636             :     !
    2637             :     !           wg    - weigths of the 30-point gauss rule
    2638             :     !
    2639             :     !
    2640             :     ! gauss quadrature weights and kronron quadrature abscissae and weights
    2641             :     ! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
    2642             :     ! bell labs, nov. 1981.
    2643             :     !
    2644             :     data wg  (  1) / 0.007968192496166605615465883474674_f /
    2645             :     data wg  (  2) / 0.018466468311090959142302131912047_f /
    2646             :     data wg  (  3) / 0.028784707883323369349719179611292_f /
    2647             :     data wg  (  4) / 0.038799192569627049596801936446348_f /
    2648             :     data wg  (  5) / 0.048402672830594052902938140422808_f /
    2649             :     data wg  (  6) / 0.057493156217619066481721689402056_f /
    2650             :     data wg  (  7) / 0.065974229882180495128128515115962_f /
    2651             :     data wg  (  8) / 0.073755974737705206268243850022191_f /
    2652             :     data wg  (  9) / 0.080755895229420215354694938460530_f /
    2653             :     data wg  ( 10) / 0.086899787201082979802387530715126_f /
    2654             :     data wg  ( 11) / 0.092122522237786128717632707087619_f /
    2655             :     data wg  ( 12) / 0.096368737174644259639468626351810_f /
    2656             :     data wg  ( 13) / 0.099593420586795267062780282103569_f /
    2657             :     data wg  ( 14) / 0.101762389748405504596428952168554_f /
    2658             :     data wg  ( 15) / 0.102852652893558840341285636705415_f /
    2659             : 
    2660             :     data xgk (  1) / 0.999484410050490637571325895705811_f /
    2661             :     data xgk (  2) / 0.996893484074649540271630050918695_f /
    2662             :     data xgk (  3) / 0.991630996870404594858628366109486_f /
    2663             :     data xgk (  4) / 0.983668123279747209970032581605663_f /
    2664             :     data xgk (  5) / 0.973116322501126268374693868423707_f /
    2665             :     data xgk (  6) / 0.960021864968307512216871025581798_f /
    2666             :     data xgk (  7) / 0.944374444748559979415831324037439_f /
    2667             :     data xgk (  8) / 0.926200047429274325879324277080474_f /
    2668             :     data xgk (  9) / 0.905573307699907798546522558925958_f /
    2669             :     data xgk ( 10) / 0.882560535792052681543116462530226_f /
    2670             :     data xgk ( 11) / 0.857205233546061098958658510658944_f /
    2671             :     data xgk ( 12) / 0.829565762382768397442898119732502_f /
    2672             :     data xgk ( 13) / 0.799727835821839083013668942322683_f /
    2673             :     data xgk ( 14) / 0.767777432104826194917977340974503_f /
    2674             :     data xgk ( 15) / 0.733790062453226804726171131369528_f /
    2675             :     data xgk ( 16) / 0.697850494793315796932292388026640_f /
    2676             :     data xgk ( 17) / 0.660061064126626961370053668149271_f /
    2677             :     data xgk ( 18) / 0.620526182989242861140477556431189_f /
    2678             :     data xgk ( 19) / 0.579345235826361691756024932172540_f /
    2679             :     data xgk ( 20) / 0.536624148142019899264169793311073_f /
    2680             :     data xgk ( 21) / 0.492480467861778574993693061207709_f /
    2681             :     data xgk ( 22) / 0.447033769538089176780609900322854_f /
    2682             :     data xgk ( 23) / 0.400401254830394392535476211542661_f /
    2683             :     data xgk ( 24) / 0.352704725530878113471037207089374_f /
    2684             :     data xgk ( 25) / 0.304073202273625077372677107199257_f /
    2685             :     data xgk ( 26) / 0.254636926167889846439805129817805_f /
    2686             :     data xgk ( 27) / 0.204525116682309891438957671002025_f /
    2687             :     data xgk ( 28) / 0.153869913608583546963794672743256_f /
    2688             :     data xgk ( 29) / 0.102806937966737030147096751318001_f /
    2689             :     data xgk ( 30) / 0.051471842555317695833025213166723_f /
    2690             :     data xgk ( 31) / 0.000000000000000000000000000000000_f /
    2691             : 
    2692             :     data wgk (  1) / 0.001389013698677007624551591226760_f /
    2693             :     data wgk (  2) / 0.003890461127099884051267201844516_f /
    2694             :     data wgk (  3) / 0.006630703915931292173319826369750_f /
    2695             :     data wgk (  4) / 0.009273279659517763428441146892024_f /
    2696             :     data wgk (  5) / 0.011823015253496341742232898853251_f /
    2697             :     data wgk (  6) / 0.014369729507045804812451432443580_f /
    2698             :     data wgk (  7) / 0.016920889189053272627572289420322_f /
    2699             :     data wgk (  8) / 0.019414141193942381173408951050128_f /
    2700             :     data wgk (  9) / 0.021828035821609192297167485738339_f /
    2701             :     data wgk ( 10) / 0.024191162078080601365686370725232_f /
    2702             :     data wgk ( 11) / 0.026509954882333101610601709335075_f /
    2703             :     data wgk ( 12) / 0.028754048765041292843978785354334_f /
    2704             :     data wgk ( 13) / 0.030907257562387762472884252943092_f /
    2705             :     data wgk ( 14) / 0.032981447057483726031814191016854_f /
    2706             :     data wgk ( 15) / 0.034979338028060024137499670731468_f /
    2707             :     data wgk ( 16) / 0.036882364651821229223911065617136_f /
    2708             :     data wgk ( 17) / 0.038678945624727592950348651532281_f /
    2709             :     data wgk ( 18) / 0.040374538951535959111995279752468_f /
    2710             :     data wgk ( 19) / 0.041969810215164246147147541285970_f /
    2711             :     data wgk ( 20) / 0.043452539701356069316831728117073_f /
    2712             :     data wgk ( 21) / 0.044814800133162663192355551616723_f /
    2713             :     data wgk ( 22) / 0.046059238271006988116271735559374_f /
    2714             :     data wgk ( 23) / 0.047185546569299153945261478181099_f /
    2715             :     data wgk ( 24) / 0.048185861757087129140779492298305_f /
    2716             :     data wgk ( 25) / 0.049055434555029778887528165367238_f /
    2717             :     data wgk ( 26) / 0.049795683427074206357811569379942_f /
    2718             :     data wgk ( 27) / 0.050405921402782346840893085653585_f /
    2719             :     data wgk ( 28) / 0.050881795898749606492297473049805_f /
    2720             :     data wgk ( 29) / 0.051221547849258772170656282604944_f /
    2721             :     data wgk ( 30) / 0.051426128537459025933862879215781_f /
    2722             :     data wgk ( 31) / 0.051494729429451567558340433647099_f /
    2723             :     
    2724             :     !           list of major variables
    2725             :     !           -----------------------
    2726             :     !
    2727             :     !           centr  - mid point of the interval
    2728             :     !           hlgth  - half-length of the interval
    2729             :     !           dabsc  - abscissa
    2730             :     !           fval*  - function value
    2731             :     !           resg   - result of the 30-point gauss rule
    2732             :     !           resk   - result of the 61-point kronrod rule
    2733             :     !           reskh  - approximation to the mean value of f
    2734             :     !                    over (a,b), i.e. to i/(b-a)
    2735             :     !
    2736             :     !           machine dependent constants
    2737             :     !           ---------------------------
    2738             :     !
    2739             :     !           epmach is the largest relative spacing.
    2740             :     ! a        uflow is the smallest positive magnitude.
    2741             :     !
    2742           0 :     call sd1mach(4,epmach,ier)
    2743           0 :     if(ier.eq.9) return
    2744           0 :     call sd1mach(1,uflow,ier)
    2745           0 :     if(ier.eq.9) return
    2746             : 
    2747             :     !epmach = d1mach(4)
    2748             :     !uflow = d1mach(1)
    2749             : 
    2750           0 :     centr = 0.5_f*(b+a)
    2751           0 :     hlgth = 0.5_f*(b-a)
    2752           0 :     dhlgth = dabs(hlgth)
    2753             :     !
    2754             :     !           compute the 61-point kronrod approximation to the
    2755             :     !           integral, and estimate the absolute error.
    2756             :     !
    2757             :     !***first executable statement  dqk61
    2758           0 :     resg = 0.0_f
    2759           0 :     fc = fx(centr, fx_vars)
    2760           0 :     resk = wgk(31)*fc
    2761           0 :     resabs = dabs(resk)
    2762           0 :     do j=1,15
    2763           0 :       jtw = j*2
    2764           0 :       dabsc = hlgth*xgk(jtw)
    2765           0 :       fval1 = fx(centr-dabsc, fx_vars)
    2766           0 :       fval2 = fx(centr+dabsc, fx_vars)
    2767           0 :       fv1(jtw) = fval1
    2768           0 :       fv2(jtw) = fval2
    2769           0 :       fsum = fval1+fval2
    2770           0 :       resg = resg+wg(j)*fsum
    2771           0 :       resk = resk+wgk(jtw)*fsum
    2772           0 :       resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
    2773             :     end do
    2774           0 :     do j=1,15
    2775           0 :       jtwm1 = j*2-1
    2776           0 :       dabsc = hlgth*xgk(jtwm1)
    2777           0 :       fval1 = fx(centr-dabsc, fx_vars)
    2778           0 :       fval2 = fx(centr+dabsc, fx_vars)
    2779           0 :       fv1(jtwm1) = fval1
    2780           0 :       fv2(jtwm1) = fval2
    2781           0 :       fsum = fval1+fval2
    2782           0 :       resk = resk+wgk(jtwm1)*fsum
    2783           0 :       resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
    2784             :     end do
    2785           0 :     reskh = resk*0.5_f
    2786           0 :     resasc = wgk(31)*dabs(fc-reskh)
    2787           0 :     do j=1,30
    2788           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    2789             :     end do
    2790           0 :     result = resk*hlgth
    2791           0 :     resabs = resabs*dhlgth
    2792           0 :     resasc = resasc*dhlgth
    2793           0 :     abserr = dabs((resk-resg)*hlgth)
    2794           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0.0_f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    2795           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    2796             : 999 return
    2797             :   end subroutine dqk61
    2798             : 
    2799             :   !!
    2800             :   !!***begin prologue  dqpsrt
    2801             :   !!***refer to  dqage,dqagie,dqagpe,dqawse
    2802             :   !!***routines called  (none)
    2803             :   !!***revision date  130319   (yymmdd)
    2804             :   !!***keywords  sequential sorting
    2805             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    2806             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    2807             :   !!***purpose  this routine maintains the descending ordering in the
    2808             :   !!            list of the local error estimated resulting from the
    2809             :   !!            interval subdivision process. at each call two error
    2810             :   !!            estimates are inserted using the sequential search
    2811             :   !!            method, top-down for the largest error estimate and
    2812             :   !!            bottom-up for the smallest error estimate.
    2813             :   !!***description
    2814             :   !!
    2815             :   !!           ordering routine
    2816             :   !!           standard fortran subroutine
    2817             :   !!           double precision version
    2818             :   !!
    2819             :   !!           parameters (meaning at output)
    2820             :   !!              limit  - integer
    2821             :   !!                       maximum number of error estimates the list
    2822             :   !!                       can contain
    2823             :   !!
    2824             :   !!              last   - integer
    2825             :   !!                       number of error estimates currently in the list
    2826             :   !!
    2827             :   !!              maxerr - integer
    2828             :   !!                       maxerr points to the nrmax-th largest error
    2829             :   !!                       estimate currently in the list
    2830             :   !!
    2831             :   !!              ermax  - double precision
    2832             :   !!                       nrmax-th largest error estimate
    2833             :   !!                       ermax = elist(maxerr)
    2834             :   !!
    2835             :   !!              elist  - double precision
    2836             :   !!                       vector of dimension last containing
    2837             :   !!                       the error estimates
    2838             :   !!
    2839             :   !!              iord   - integer
    2840             :   !!                       vector of dimension last, the first k elements
    2841             :   !!                       of which contain pointers to the error
    2842             :   !!                       estimates, such that
    2843             :   !!                       elist(iord(1)),...,  elist(iord(k))
    2844             :   !!                       form a decreasing sequence, with
    2845             :   !!                       k = last if last.le.(limit/2+2), and
    2846             :   !!                       k = limit+1-last otherwise
    2847             :   !!
    2848             :   !!              nrmax  - integer
    2849             :   !!                       maxerr = iord(nrmax)
    2850             :   !!
    2851             :   !!***end prologue  dqpsrt
    2852             :   !!
    2853           0 :   subroutine dqpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
    2854             : 
    2855             :     ! Arguments
    2856             :     integer :: limit
    2857             :     integer :: last
    2858             :     integer :: maxerr
    2859             :     real(kind=f) :: ermax
    2860             :     real(kind=f) :: elist(last)
    2861             :     integer :: iord(last)
    2862             :     integer :: nrmax
    2863             : 
    2864             :     ! Local declarations
    2865             :     real(kind=f) :: errmax, errmin
    2866             :     integer :: i, ibeg, ido, isucc, j, jbnd, jupbn, k
    2867             : 
    2868             :     !
    2869             :     !           check whether the list contains more than
    2870             :     !           two error estimates.
    2871             :     !
    2872             :     !***first executable statement  dqpsrt
    2873           0 :     if(last.gt.2) go to 10
    2874           0 :     iord(1) = 1
    2875           0 :     iord(2) = 2
    2876           0 :     go to 90
    2877             :     !
    2878             :     !           this part of the routine is only executed if, due to a
    2879             :     !           difficult integrand, subdivision increased the error
    2880             :     !           estimate. in the normal case the insert procedure should
    2881             :     !           start after the nrmax-th largest error estimate.
    2882             :     !
    2883           0 :  10 errmax = elist(maxerr)
    2884           0 :     if(nrmax.eq.1) go to 30
    2885             :     ido = nrmax-1
    2886           0 :     do i = 1,ido
    2887           0 :       isucc = iord(nrmax-1)
    2888             :       ! ***jump out of do-loop
    2889           0 :       if(errmax.le.elist(isucc)) go to 30
    2890           0 :       iord(nrmax) = isucc
    2891           0 :       nrmax = nrmax-1
    2892             :     end do
    2893             :     !
    2894             :     !           compute the number of elements in the list to be maintained
    2895             :     !           in descending order. this number depends on the number of
    2896             :     !           subdivisions still allowed.
    2897             :     !
    2898           0 :  30 jupbn = last
    2899           0 :     if(last.gt.(limit/2+2)) jupbn = limit+3-last
    2900           0 :     errmin = elist(last)
    2901             :     !
    2902             :     !           insert errmax by traversing the list top-down,
    2903             :     !           starting comparison from the element elist(iord(nrmax+1)).
    2904             :     !
    2905           0 :     jbnd = jupbn-1
    2906           0 :     ibeg = nrmax+1
    2907           0 :     if(ibeg.gt.jbnd) go to 50
    2908           0 :     do i=ibeg,jbnd
    2909           0 :       isucc = iord(i)
    2910             :       ! ***jump out of do-loop
    2911           0 :       if(errmax.ge.elist(isucc)) go to 60
    2912           0 :       iord(i-1) = isucc
    2913             :     end do
    2914           0 :  50 iord(jbnd) = maxerr
    2915           0 :     iord(jupbn) = last
    2916           0 :     go to 90
    2917             :     !
    2918             :     !           insert errmin by traversing the list bottom-up.
    2919             :     !
    2920           0 :  60 iord(i-1) = maxerr
    2921           0 :     k = jbnd
    2922           0 :     do j=i,jbnd
    2923           0 :       isucc = iord(k)
    2924             :       ! ***jump out of do-loop
    2925           0 :       if(errmin.lt.elist(isucc)) go to 80
    2926           0 :       iord(k+1) = isucc
    2927           0 :       k = k-1
    2928             :     end do
    2929           0 :     iord(i) = last
    2930           0 :     go to 90
    2931           0 :  80 iord(k+1) = last
    2932             :     !
    2933             :     !           set maxerr and ermax.
    2934             :     !
    2935           0 :  90 maxerr = iord(nrmax)
    2936           0 :     ermax = elist(maxerr)
    2937           0 :     return
    2938             :   end subroutine dqpsrt
    2939             : 
    2940             :   !!
    2941             :   !!***begin prologue  dqk15i
    2942             :   !!***date written   800101   (yymmdd)
    2943             :   !!***revision date  130319   (yymmdd)
    2944             :   !!***category no.  h2a3a2,h2a4a2
    2945             :   !!***keywords  15-point transformed gauss-kronrod rules
    2946             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    2947             :   !!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
    2948             :   !!***purpose  the original (infinite integration range is mapped
    2949             :   !!            onto the interval (0,1) and (a,b) is a part of (0,1).
    2950             :   !!            it is the purpose to compute
    2951             :   !!            i = integral of transformed integrand over (a,b),
    2952             :   !!            j = integral of abs(transformed integrand) over (a,b).
    2953             :   !!***description
    2954             :   !!
    2955             :   !!           integration rule
    2956             :   !!           standard fortran subroutine
    2957             :   !!           double precision version
    2958             :   !!
    2959             :   !!           parameters
    2960             :   !!            on entry
    2961             :   !!              fx     - double precision
    2962             :   !!                       fuction subprogram defining the integrand
    2963             :   !!                       function f(x). the actual name for f needs to be
    2964             :   !!                       declared e x t e r n a l in the calling program.
    2965             :   !!
    2966             :   !!              fx_vars- structure containing variables need for integration
    2967             :   !!                       specific to fractal meanfield scattering code!
    2968             :   !!
    2969             :   !!              boun   - double precision
    2970             :   !!                       finite bound of original integration
    2971             :   !!                       range (set to zero if inf = +2)
    2972             :   !!
    2973             :   !!              inf    - integer
    2974             :   !!                       if inf = -1, the original interval is
    2975             :   !!                                   (-infinity,bound),
    2976             :   !!                       if inf = +1, the original interval is
    2977             :   !!                                   (bound,+infinity),
    2978             :   !!                       if inf = +2, the original interval is
    2979             :   !!                                   (-infinity,+infinity) and
    2980             :   !!                       the integral is computed as the sum of two
    2981             :   !!                       integrals, one over (-infinity,0) and one over
    2982             :   !!                       (0,+infinity).
    2983             :   !!
    2984             :   !!              a      - double precision
    2985             :   !!                       lower limit for integration over subrange
    2986             :   !!                       of (0,1)
    2987             :   !!
    2988             :   !!              b      - double precision
    2989             :   !!                       upper limit for integration over subrange
    2990             :   !!                       of (0,1)
    2991             :   !!
    2992             :   !!            on return
    2993             :   !!              result - double precision
    2994             :   !!                       approximation to the integral i
    2995             :   !!                       result is computed by applying the 15-point
    2996             :   !!                       kronrod rule(resk) obtained by optimal addition
    2997             :   !!                       of abscissae to the 7-point gauss rule(resg).
    2998             :   !!
    2999             :   !!              abserr - double precision
    3000             :   !!                       estimate of the modulus of the absolute error,
    3001             :   !!                       which should equal or exceed abs(i-result)
    3002             :   !!
    3003             :   !!              resabs - double precision
    3004             :   !!                       approximation to the integral j
    3005             :   !!
    3006             :   !!              resasc - double precision
    3007             :   !!                       approximation to the integral of
    3008             :   !!                       abs((transformed integrand)-i/(b-a)) over (a,b)
    3009             :   !!
    3010             :   !!              ier    - integer
    3011             :   !!                       ier = 0 normal and reliable termination of the
    3012             :   !!                               routine. it is assumed that the requested
    3013             :   !!                               accuracy has been achieved.
    3014             :   !!                       ier.gt.0 abnormal termination of the routine. the
    3015             :   !!                                estimates for result and error are less
    3016             :   !!                                reliable. it is assumed that the requested
    3017             :   !!                                accuracy has not been achieved.
    3018             :   !!
    3019             :   !!***references  (none)
    3020             :   !!***routines called  sd1mach
    3021             :   !!***end prologue  dqk15i
    3022           0 :   subroutine dqk15i(fx,fx_vars,boun,inf,a,b,result,abserr,resabs,resasc, ier)
    3023             : 
    3024             :     ! Arguments
    3025             :     interface
    3026             :       function fx(centr, vars)
    3027             :         use carma_precision_mod, only : f
    3028             :         use adgaquad_types_mod
    3029             :         real(kind=f), intent(in) :: centr
    3030             :         type(adgaquad_vars_type), intent(inout) :: vars
    3031             :         real(kind=f)             :: fx
    3032             :       end function fx
    3033             :     end interface
    3034             :     type(adgaquad_vars_type) :: fx_vars
    3035             :     real(kind=f) :: boun
    3036             :     integer :: inf
    3037             :     real(kind=f) :: a 
    3038             :     real(kind=f) :: b
    3039             :     real(kind=f) :: result
    3040             :     real(kind=f) :: abserr
    3041             :     real(kind=f) :: resabs
    3042             :     real(kind=f) :: resasc
    3043             :     integer :: ier
    3044             : 
    3045             :     ! Local declarations
    3046             :     real(kind=f) :: absc, absc1, absc2, centr, dabs, dinf
    3047             :     real(kind=f) :: dmax1, dmin1, epmach, fc, fsum, fval1, fval2, fv1(7) ,fv2(7), hlgth
    3048             :     real(kind=f) :: resg, resk, reskh, tabsc1, tabsc2, uflow, wg(8), wgk(8), xgk(8)
    3049             :     integer :: j
    3050             : 
    3051             :     !
    3052             :     !           the abscissae and weights are supplied for the interval
    3053             :     !           (-1,1).  because of symmetry only the positive abscissae and
    3054             :     !           their corresponding weights are given.
    3055             :     !
    3056             :     !           xgk    - abscissae of the 15-point kronrod rule
    3057             :     !                    xgk(2), xgk(4), ... abscissae of the 7-point
    3058             :     !                    gauss rule
    3059             :     !                    xgk(1), xgk(3), ...  abscissae which are optimally
    3060             :     !                    added to the 7-point gauss rule
    3061             :     !
    3062             :     !           wgk    - weights of the 15-point kronrod rule
    3063             :     !
    3064             :     !           wg     - weights of the 7-point gauss rule, corresponding
    3065             :     !                    to the abscissae xgk(2), xgk(4), ...
    3066             :     !                    wg(1), wg(3), ... are set to zero.
    3067             :     !
    3068             :     data wg(1) / 0.0_f /
    3069             :     data wg(2) / 0.129484966168869693270611432679082_f /
    3070             :     data wg(3) / 0.0_f /
    3071             :     data wg(4) / 0.279705391489276667901467771423780_f /
    3072             :     data wg(5) / 0.0_f /
    3073             :     data wg(6) / 0.381830050505118944950369775488975_f /
    3074             :     data wg(7) / 0.0_f /
    3075             :     data wg(8) / 0.417959183673469387755102040816327_f /
    3076             : 
    3077             :     data xgk(1) / 0.991455371120812639206854697526329_f /
    3078             :     data xgk(2) / 0.949107912342758524526189684047851_f /
    3079             :     data xgk(3) / 0.864864423359769072789712788640926_f /
    3080             :     data xgk(4) / 0.741531185599394439863864773280788_f /
    3081             :     data xgk(5) / 0.586087235467691130294144838258730_f /
    3082             :     data xgk(6) / 0.405845151377397166906606412076961_f /
    3083             :     data xgk(7) / 0.207784955007898467600689403773245_f /
    3084             :     data xgk(8) / 0.000000000000000000000000000000000_f /
    3085             : 
    3086             :     data wgk(1) / 0.022935322010529224963732008058970_f /
    3087             :     data wgk(2) / 0.063092092629978553290700663189204_f /
    3088             :     data wgk(3) / 0.104790010322250183839876322541518_f /
    3089             :     data wgk(4) / 0.140653259715525918745189590510238_f /
    3090             :     data wgk(5) / 0.169004726639267902826583426598550_f /
    3091             :     data wgk(6) / 0.190350578064785409913256402421014_f /
    3092             :     data wgk(7) / 0.204432940075298892414161999234649_f /
    3093             :     data wgk(8) / 0.209482141084727828012999174891714_f /
    3094             :     !
    3095             :     !
    3096             :     !           list of major variables
    3097             :     !           -----------------------
    3098             :     !
    3099             :     !           centr  - mid point of the interval
    3100             :     !           hlgth  - half-length of the interval
    3101             :     !           absc*  - abscissa
    3102             :     !           tabsc* - transformed abscissa
    3103             :     !           fval*  - function value
    3104             :     !           resg   - result of the 7-point gauss formula
    3105             :     !           resk   - result of the 15-point kronrod formula
    3106             :     !           reskh  - approximation to the mean value of the transformed
    3107             :     !                    integrand over (a,b), i.e. to i/(b-a)
    3108             :     !
    3109             :     !           machine dependent constants
    3110             :     !           ---------------------------
    3111             :     !
    3112             :     !           epmach is the largest relative spacing.
    3113             :     !           uflow is the smallest positive magnitude.
    3114             :     !
    3115             :     !*** first executable statement  dqk15i
    3116           0 :     call sd1mach(4,epmach,ier)
    3117           0 :     if(ier.eq.9) return
    3118           0 :     call sd1mach(1,uflow,ier)
    3119           0 :     if(ier.eq.9) return
    3120             : 
    3121             :     !epmach = d1mach(4)
    3122             :     !uflow = d1mach(1)
    3123           0 :     dinf = min0(1,inf)
    3124             : 
    3125           0 :     centr = 0.5_f*(a+b)
    3126           0 :     hlgth = 0.5_f*(b-a)
    3127           0 :     tabsc1 = boun+dinf*(0.1e1_f-centr)/centr
    3128           0 :     fval1 = fx(tabsc1, fx_vars)
    3129           0 :     if(inf.eq.2) fval1 = fval1+fx(-tabsc1, fx_vars)
    3130           0 :     fc = (fval1/centr)/centr
    3131             :     !
    3132             :     !           compute the 15-point kronrod approximation to
    3133             :     !           the integral, and estimate the error.
    3134             :     !
    3135           0 :     resg = wg(8)*fc
    3136           0 :     resk = wgk(8)*fc
    3137           0 :     resabs = dabs(resk)
    3138           0 :     do j=1,7
    3139           0 :       absc = hlgth*xgk(j)
    3140           0 :       absc1 = centr-absc
    3141           0 :       absc2 = centr+absc
    3142           0 :       tabsc1 = boun+dinf*(1.0_f-absc1)/absc1
    3143           0 :       tabsc2 = boun+dinf*(1.0_f-absc2)/absc2
    3144           0 :       fval1 = fx(tabsc1, fx_vars)
    3145           0 :       fval2 = fx(tabsc2, fx_vars)
    3146           0 :       if(inf.eq.2) fval1 = fval1+fx(-tabsc1, fx_vars)
    3147           0 :       if(inf.eq.2) fval2 = fval2+fx(-tabsc2, fx_vars)
    3148           0 :       fval1 = (fval1/absc1)/absc1
    3149           0 :       fval2 = (fval2/absc2)/absc2
    3150           0 :       fv1(j) = fval1
    3151           0 :       fv2(j) = fval2
    3152           0 :       fsum = fval1+fval2
    3153           0 :       resg = resg+wg(j)*fsum
    3154           0 :       resk = resk+wgk(j)*fsum
    3155           0 :       resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2))
    3156             :     end do
    3157           0 :     reskh = resk*0.5_f
    3158           0 :     resasc = wgk(8)*dabs(fc-reskh)
    3159           0 :     do j=1,7
    3160           0 :       resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
    3161             :     end do
    3162           0 :     result = resk*hlgth
    3163           0 :     resasc = resasc*hlgth
    3164           0 :     resabs = resabs*hlgth
    3165           0 :     abserr = dabs((resk-resg)*hlgth)
    3166           0 :     if(resasc.ne.0.0_f.and.abserr.ne.0._f) abserr = resasc*dmin1(0.1e1_f,(0.2e3_f*abserr/resasc)**1.5_f)
    3167           0 :     if(resabs.gt.uflow/(0.5e2_f*epmach)) abserr = dmax1((epmach*0.5e2_f)*resabs,abserr)
    3168             : 999 return
    3169             :   end subroutine dqk15i
    3170             : 
    3171             :   !!
    3172             :   !!***begin prologue  dqelg
    3173             :   !!***refer to  dqagie,dqagoe,dqagpe,dqagse
    3174             :   !!***routines called  sd1mach
    3175             :   !!***revision date  130319   (yymmdd)
    3176             :   !!***keywords  epsilon algorithm, convergence acceleration,
    3177             :   !!             extrapolation
    3178             :   !!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
    3179             :   !!           de doncker,elise,appl. math & progr. div. - k.u.leuven
    3180             :   !!***purpose  the routine determines the limit of a given sequence of
    3181             :   !!            approximations, by means of the epsilon algorithm of
    3182             :   !!            p.wynn. an estimate of the absolute error is also given.
    3183             :   !!            the condensed epsilon table is computed. only those
    3184             :   !!            elements needed for the computation of the next diagonal
    3185             :   !!            are preserved.
    3186             :   !!***description
    3187             :   !!
    3188             :   !!           epsilon algorithm
    3189             :   !!           standard fortran subroutine
    3190             :   !!           double precision version
    3191             :   !!
    3192             :   !!           parameters
    3193             :   !!              n      - integer
    3194             :   !!                       epstab(n) contains the new element in the
    3195             :   !!                       first column of the epsilon table.
    3196             :   !!
    3197             :   !!              epstab - double precision
    3198             :   !!                       vector of dimension 52 containing the elements
    3199             :   !!                       of the two lower diagonals of the triangular
    3200             :   !!                       epsilon table. the elements are numbered
    3201             :   !!                       starting at the right-hand corner of the
    3202             :   !!                       triangle.
    3203             :   !!
    3204             :   !!              result - double precision
    3205             :   !!                       resulting approximation to the integral
    3206             :   !!
    3207             :   !!              abserr - double precision
    3208             :   !!                       estimate of the absolute error computed from
    3209             :   !!                       result and the 3 previous results
    3210             :   !!
    3211             :   !!              res3la - double precision
    3212             :   !!                       vector of dimension 3 containing the last 3
    3213             :   !!                       results
    3214             :   !!
    3215             :   !!              nres   - integer
    3216             :   !!                       number of calls to the routine
    3217             :   !!                       (should be zero at first call)
    3218             :   !!
    3219             :   !!              ier    - integer
    3220             :   !!                       ier = 0 normal and reliable termination of the
    3221             :   !!                               routine. it is assumed that the requested
    3222             :   !!                               accuracy has been achieved.
    3223             :   !!                       ier.gt.0 abnormal termination of the routine. the
    3224             :   !!                                estimates for result and error are less
    3225             :   !!                                reliable. it is assumed that the requested
    3226             :   !!                                accuracy has not been achieved.
    3227             :   !!
    3228             :   !!***end prologue  dqelg
    3229             :   !!
    3230           0 :   subroutine dqelg(n,epstab,result,abserr,res3la,nres,ier)
    3231             : 
    3232             :     ! Arguments
    3233             :     integer :: n
    3234             :     real(kind=f) :: epstab(52)
    3235             :     real(kind=f) :: result
    3236             :     real(kind=f) :: abserr
    3237             :     real(kind=f) :: res3la(3)
    3238             :     integer :: nres  
    3239             :     integer :: ier
    3240             : 
    3241             :     ! Local declarations
    3242             :     real(kind=f) :: dabs, delta1, delta2, delta3, dmax1
    3243             :     real(kind=f) :: epmach, epsinf, error, err1, err2, err3, e0, e1, e1abs, e2, e3
    3244             :     real(kind=f) :: oflow, res, ss, tol1, tol2, tol3
    3245             :     integer :: i, ib, ib2, ie, indx, k1, k2, k3, limexp, newelm, num
    3246             :     !
    3247             :     !           list of major variables
    3248             :     !           -----------------------
    3249             :     !
    3250             :     !           e0     - the 4 elements on which the computation of a new
    3251             :     !           e1       element in the epsilon table is based
    3252             :     !           e2
    3253             :     !           e3                 e0
    3254             :     !                        e3    e1    new
    3255             :     !                              e2
    3256             :     !           newelm - number of elements to be computed in the new
    3257             :     !                    diagonal
    3258             :     !           error  - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
    3259             :     !           result - the element in the new diagonal with least value
    3260             :     !                    of error
    3261             :     !
    3262             :     !           machine dependent constants
    3263             :     !           ---------------------------
    3264             :     !
    3265             :     !           epmach is the largest relative spacing.
    3266             :     !           oflow is the largest positive magnitude.
    3267             :     !           limexp is the maximum number of elements the epsilon
    3268             :     !           table can contain. if this number is reached, the upper
    3269             :     !           diagonal of the epsilon table is deleted.
    3270             :     !
    3271             :     !***first executable statement  dqelg
    3272           0 :     call sd1mach(4,epmach,ier)
    3273           0 :     if(ier.eq.9) return
    3274           0 :     call sd1mach(2,oflow,ier)
    3275           0 :     if(ier.eq.9) return
    3276             : 
    3277             :     !epmach = d1mach(4)
    3278             :     !oflow = d1mach(2)
    3279           0 :     nres = nres+1
    3280           0 :     abserr = oflow
    3281           0 :     result = epstab(n)
    3282           0 :     if(n.lt.3) go to 100
    3283           0 :     limexp = 50
    3284           0 :     epstab(n+2) = epstab(n)
    3285           0 :     newelm = (n-1)/2
    3286           0 :     epstab(n) = oflow
    3287           0 :     num = n
    3288           0 :     k1 = n
    3289           0 :     do 40 i = 1,newelm
    3290           0 :       k2 = k1-1
    3291           0 :       k3 = k1-2
    3292           0 :       res = epstab(k1+2)
    3293           0 :       e0 = epstab(k3)
    3294           0 :       e1 = epstab(k2)
    3295           0 :       e2 = res
    3296           0 :       e1abs = dabs(e1)
    3297           0 :       delta2 = e2-e1
    3298           0 :       err2 = dabs(delta2)
    3299           0 :       tol2 = dmax1(dabs(e2),e1abs)*epmach
    3300           0 :       delta3 = e1-e0
    3301           0 :       err3 = dabs(delta3)
    3302           0 :       tol3 = dmax1(e1abs,dabs(e0))*epmach
    3303           0 :       if(err2.gt.tol2.or.err3.gt.tol3) go to 10
    3304             :       !
    3305             :       !           if e0, e1 and e2 are equal to within machine
    3306             :       !           accuracy, convergence is assumed.
    3307             :       !           result = e2
    3308             :       !           abserr = abs(e1-e0)+abs(e2-e1)
    3309             :       !
    3310           0 :       result = res
    3311           0 :       abserr = err2+err3
    3312             :       ! ***jump out of do-loop
    3313           0 :       go to 100
    3314           0 :  10   e3 = epstab(k1)
    3315           0 :       epstab(k1) = e1
    3316           0 :       delta1 = e1-e3
    3317           0 :       err1 = dabs(delta1)
    3318           0 :       tol1 = dmax1(e1abs,dabs(e3))*epmach
    3319             :       !
    3320             :       !           if two elements are very close to each other, omit
    3321             :       !           a part of the table by adjusting the value of n
    3322             :       !
    3323           0 :       if(err1.le.tol1.or.err2.le.tol2.or.err3.le.tol3) go to 20
    3324           0 :       ss = 0.1e1_f/delta1+0.1e1_f/delta2-0.1e1_f/delta3
    3325           0 :       epsinf = dabs(ss*e1)
    3326             :       !
    3327             :       !           test to detect irregular behaviour in the table, and
    3328             :       !           eventually omit a part of the table adjusting the value
    3329             :       !           of n.
    3330             :       !
    3331           0 :       if(epsinf.gt.0.1e-3_f) go to 30
    3332           0 :  20   n = i+i-1
    3333             :       ! ***jump out of do-loop
    3334           0 :       go to 50
    3335             :       !
    3336             :       !           compute a new element and eventually adjust
    3337             :       !           the value of result.
    3338             :       !
    3339           0 :  30   res = e1+0.1e1_f/ss
    3340           0 :       epstab(k1) = res
    3341           0 :       k1 = k1-2
    3342           0 :       error = err2+dabs(res-e2)+err3
    3343           0 :       if(error.gt.abserr) go to 40
    3344           0 :       abserr = error
    3345           0 :       result = res
    3346           0 :  40 continue
    3347             :     !
    3348             :     !           shift the table.
    3349             :     !
    3350           0 :  50 if(n.eq.limexp) n = 2*(limexp/2)-1
    3351           0 :     ib = 1
    3352           0 :     if((num/2)*2.eq.num) ib = 2
    3353           0 :     ie = newelm+1
    3354           0 :     do 60 i=1,ie
    3355           0 :       ib2 = ib+2
    3356           0 :       epstab(ib) = epstab(ib2)
    3357           0 :       ib = ib2
    3358           0 :  60 continue
    3359           0 :     if(num.eq.n) go to 80
    3360           0 :     indx = num-n+1
    3361           0 :     do 70 i = 1,n
    3362           0 :       epstab(i)= epstab(indx)
    3363           0 :       indx = indx+1
    3364           0 :  70 continue
    3365           0 :  80 if(nres.ge.4) go to 90
    3366           0 :     res3la(nres) = result
    3367           0 :     abserr = oflow
    3368           0 :     go to 100
    3369             :     !
    3370             :     !           compute error estimate
    3371             :     !
    3372           0 :  90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))+dabs(result-res3la(1))
    3373           0 :     res3la(1) = res3la(2)
    3374           0 :     res3la(2) = res3la(3)
    3375           0 :     res3la(3) = result
    3376           0 : 100 abserr = dmax1(abserr,0.5e1_f*epmach*dabs(result))
    3377           0 : 999 return
    3378             :   end subroutine dqelg
    3379             :   
    3380             :   !!
    3381             :   !! *********************************************************
    3382             :   !! taken from BLAS library
    3383             :   !!  (http://netlib.bell-labs.com/netlib/blas)
    3384             :   !! *********************************************************
    3385           0 :   SUBROUTINE SD1MACH(I,D1MACH_OUT,IER)
    3386             :     INTEGER, INTENT(in) :: I
    3387             :     REAL(kind=f), INTENT(out) :: D1MACH_OUT
    3388             :     INTEGER, INTENT(out) :: IER
    3389             :   !
    3390             :   !  DOUBLE-PRECISION MACHINE CONSTANTS
    3391             :   !  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
    3392             :   !  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
    3393             :   !  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
    3394             :   !  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
    3395             :   !  D1MACH( 5) = LOG10(B)
    3396             :   !
    3397             :     INTEGER :: SMALL(2)
    3398             :     INTEGER :: LARGE(2)
    3399             :     INTEGER :: RIGHT(2)
    3400             :     INTEGER :: DIVER(2)
    3401             :     INTEGER :: LOG10(2)
    3402             :     INTEGER :: SC, CRAY1(38), J
    3403             :     SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC
    3404             :     REAL(kind=f) :: DMACH(5)
    3405             :     EQUIVALENCE (DMACH(1),SMALL(1))
    3406             :     EQUIVALENCE (DMACH(2),LARGE(1))
    3407             :     EQUIVALENCE (DMACH(3),RIGHT(1))
    3408             :     EQUIVALENCE (DMACH(4),DIVER(1))
    3409             :     EQUIVALENCE (DMACH(5),LOG10(1))
    3410             :     !  THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES.
    3411             :     !  R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF
    3412             :     !  D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR
    3413             :     !  MANY MACHINES YET.
    3414             :     !  TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1
    3415             :     !  ON THE NEXT LINE
    3416             :     DATA SC/0/
    3417             :     !  AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW.
    3418             :     !  CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY
    3419             :     !          mail netlib@research.bell-labs.com
    3420             :     !          send old1mach from blas
    3421             :     !  PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com.
    3422             :     !
    3423             :     !     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
    3424             :     !      DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
    3425             :     !      DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
    3426             :     !      DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
    3427             :     !      DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
    3428             :     !      DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/
    3429             :     !
    3430             :     !     MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING
    3431             :     !     32-BIT INTEGERS.
    3432             :     !      DATA SMALL(1),SMALL(2) /    8388608,           0 /
    3433             :     !      DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
    3434             :     !      DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
    3435             :     !      DATA DIVER(1),DIVER(2) /  620756992,           0 /
    3436             :     !      DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/
    3437             :     !
    3438             :     !  MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
    3439             :     !      DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
    3440             :     !      DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
    3441             :     !      DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
    3442             :     !      DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
    3443             :     !      DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/
    3444             :     !
    3445             :     !     ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES.
    3446           0 :     IER = 0
    3447           0 :     IF (SC .NE. 987) THEN
    3448           0 :       DMACH(1) = 1.e13_f
    3449             :       IF (      SMALL(1) .EQ. 1117925532 .AND. SMALL(2) .EQ. -448790528) THEN
    3450             :       !           *** IEEE BIG ENDIAN ***
    3451             :         SMALL(1) = 1048576
    3452             :         SMALL(2) = 0
    3453             :         LARGE(1) = 2146435071
    3454             :         LARGE(2) = -1
    3455             :         RIGHT(1) = 1017118720
    3456             :         RIGHT(2) = 0
    3457             :         DIVER(1) = 1018167296
    3458             :         DIVER(2) = 0
    3459             :         LOG10(1) = 1070810131
    3460             :         LOG10(2) = 1352628735
    3461             :       ELSE IF ( SMALL(2) .EQ. 1117925532 .AND. SMALL(1) .EQ. -448790528) THEN
    3462             :       !           *** IEEE LITTLE ENDIAN ***
    3463           0 :         SMALL(2) = 1048576
    3464           0 :         SMALL(1) = 0
    3465           0 :         LARGE(2) = 2146435071
    3466           0 :         LARGE(1) = -1
    3467           0 :         RIGHT(2) = 1017118720
    3468           0 :         RIGHT(1) = 0
    3469           0 :         DIVER(2) = 1018167296
    3470           0 :         DIVER(1) = 0
    3471           0 :         LOG10(2) = 1070810131
    3472           0 :         LOG10(1) = 1352628735
    3473             :       ELSE IF ( SMALL(1) .EQ. -2065213935 .AND. SMALL(2) .EQ. 10752) THEN
    3474             :       !               *** VAX WITH D_FLOATING ***
    3475             :         SMALL(1) = 128
    3476             :         SMALL(2) = 0
    3477             :         LARGE(1) = -32769
    3478             :         LARGE(2) = -1
    3479             :         RIGHT(1) = 9344
    3480             :         RIGHT(2) = 0
    3481             :         DIVER(1) = 9472
    3482             :         DIVER(2) = 0
    3483             :         LOG10(1) = 546979738
    3484             :         LOG10(2) = -805796613
    3485             :       ELSE IF ( SMALL(1) .EQ. 1267827943 .AND. SMALL(2) .EQ. 704643072) THEN
    3486             :       !               *** IBM MAINFRAME ***
    3487             :         SMALL(1) = 1048576
    3488             :         SMALL(2) = 0
    3489             :         LARGE(1) = 2147483647
    3490             :         LARGE(2) = -1
    3491             :         RIGHT(1) = 856686592
    3492             :         RIGHT(2) = 0
    3493             :         DIVER(1) = 873463808
    3494             :         DIVER(2) = 0
    3495             :         LOG10(1) = 1091781651
    3496             :         LOG10(2) = 1352628735
    3497             :       ELSE IF ( SMALL(1) .EQ. 1120022684 .AND. SMALL(2) .EQ. -448790528) THEN
    3498             :       !           *** CONVEX C-1 ***
    3499             :         SMALL(1) = 1048576
    3500             :         SMALL(2) = 0
    3501             :         LARGE(1) = 2147483647
    3502             :         LARGE(2) = -1
    3503             :         RIGHT(1) = 1019215872
    3504             :         RIGHT(2) = 0
    3505             :         DIVER(1) = 1020264448
    3506             :         DIVER(2) = 0
    3507             :         LOG10(1) = 1072907283
    3508             :         LOG10(2) = 1352628735
    3509             :       ELSE IF ( SMALL(1) .EQ. 815547074 .AND. SMALL(2) .EQ. 58688) THEN
    3510             :       !           *** VAX G-FLOATING ***
    3511             :         SMALL(1) = 16
    3512             :         SMALL(2) = 0
    3513             :         LARGE(1) = -32769
    3514             :         LARGE(2) = -1
    3515             :         RIGHT(1) = 15552
    3516             :         RIGHT(2) = 0
    3517             :         DIVER(1) = 15568
    3518             :         DIVER(2) = 0
    3519             :         LOG10(1) = 1142112243
    3520             :         LOG10(2) = 2046775455
    3521             :       ELSE
    3522             :         DMACH(2) = 1.e27_f + 1
    3523             :         DMACH(3) = 1.e27_f
    3524             :         LARGE(2) = LARGE(2) - RIGHT(2)
    3525             :         IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN
    3526             :           CRAY1(1) = 67291416
    3527             :           DO J = 1, 20
    3528             :             CRAY1(J+1) = CRAY1(J) + CRAY1(J)
    3529             :           END DO
    3530             :           CRAY1(22) = CRAY1(21) + 321322
    3531             :           DO J = 22, 37
    3532             :             CRAY1(J+1) = CRAY1(J) + CRAY1(J)
    3533             :           END DO
    3534             :           IF (CRAY1(38) .EQ. SMALL(1)) THEN
    3535             :           !                  *** CRAY ***
    3536             :             CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0)
    3537             :             SMALL(2) = 0
    3538             :             CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215)
    3539             :             CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214)
    3540             :             CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0)
    3541             :             RIGHT(2) = 0
    3542             :             CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0)
    3543             :             DIVER(2) = 0
    3544             :             CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215)
    3545             :             CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388)
    3546             :           ELSE
    3547             :             IER=9
    3548             :           END IF
    3549             :           ELSE
    3550             :             IER=9
    3551             :           END IF
    3552             :         END IF
    3553           0 :         SC = 987
    3554             :       END IF
    3555             :       !    SANITY CHECK
    3556           0 :       IF (DMACH(4) .GE. 1.0D0) IER=9
    3557           0 :       IF (I .LT. 1 .OR. I .GT. 5) THEN
    3558           0 :         IER=9
    3559             :       END IF
    3560           0 :       D1MACH_OUT = DMACH(I)
    3561           0 :       RETURN
    3562             :   END SUBROUTINE SD1MACH
    3563             : 
    3564             :   SUBROUTINE I1MCRY(A, A1, B, C, D)
    3565             :     !*** SPECIAL COMPUTATION FOR OLD CRAY MACHINES ****
    3566             :     INTEGER A, A1, B, C, D
    3567             :     A1 = 16777216*B + C
    3568             :     A = 16777216*A1 + D
    3569             :   END SUBROUTINE I1MCRY
    3570             : 
    3571             : 
    3572             : 
    3573             : end module adgaquad_mod