!   fvn comment :
!   Modified version of the dqage quadpack routine from http://www.netlib.org/quadpack
!
!   + The external 'f' function is a 2 parameters function f(x,y). The routine
!   takes two more parameters 'g' and 'h' which are two external functions :
!   g represent the lower bound of the integral for y parameter
!   h represent the higher bound of the integral for y parameter
!   The routine compute the double integral of function f with x between a and b
!   and y between g(x) and h(x)
      subroutine dqage_2d_outer(f,a,b,g,h,epsabs,epsrel,key, &
        limit,result,abserr,neval,ier,alist,blist,rlist, &
        elist,iord,last)
!***begin prologue  dqage
!***date written   800101   (yymmdd)
!***revision date  830518   (yymmdd)
!***category no.  h2a1a1
!***keywords  automatic integrator, general-purpose,
!             integrand examinator, globally adaptive,
!             gauss-kronrod
!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
!***purpose  the routine calculates an approximation result to a given
!            definite integral   i = integral of f over (a,b),
!            hopefully satisfying following claim for accuracy
!            abs(i-reslt).le.max(epsabs,epsrel*abs(i)).
!***description
!
!        computation of a definite integral
!        standard fortran subroutine
!        double precision version
!
!        parameters
!         on entry
!            f      - double precision
!                     function subprogram defining the integrand
!                     function f(x). the actual name for f needs to be
!                     declared e x t e r n a l in the driver program.
!
!            a      - double precision
!                     lower limit of integration
!
!            b      - double precision
!                     upper limit of integration
!
!            epsabs - double precision
!                     absolute accuracy requested
!            epsrel - double precision
!                     relative accuracy requested
!                     if  epsabs.le.0
!                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
!                     the routine will end with ier = 6.
!
!            key    - integer
!                     key for choice of local integration rule
!                     a gauss-kronrod pair is used with
!                          7 - 15 points if key.lt.2,
!                         10 - 21 points if key = 2,
!                         15 - 31 points if key = 3,
!                         20 - 41 points if key = 4,
!                         25 - 51 points if key = 5,
!                         30 - 61 points if key.gt.5.
!
!            limit  - integer
!                     gives an upperbound on the number of subintervals
!                     in the partition of (a,b), limit.ge.1.
!
!         on return
!            result - double precision
!                     approximation to the integral
!
!            abserr - double precision
!                     estimate of the modulus of the absolute error,
!                     which should equal or exceed abs(i-result)
!
!            neval  - integer
!                     number of integrand evaluations
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier.gt.0 abnormal termination of the routine
!                             the estimates for result and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!            error messages
!                     ier = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more
!                             subdivisions by increasing the value
!                             of limit.
!                             however, if this yields no improvement it
!                             is rather advised to analyze the integrand
!                             in order to determine the integration
!                             difficulties. if the position of a local
!                             difficulty can be determined(e.g.
!                             singularity, discontinuity within the
!                             interval) one will probably gain from
!                             splitting up the interval at this point
!                             and calling the integrator on the
!                             subranges. if possible, an appropriate
!                             special-purpose integrator should be used
!                             which is designed for handling the type of
!                             difficulty involved.
!                         = 2 the occurrence of roundoff error is
!                             detected, which prevents the requested
!                             tolerance from being achieved.
!                         = 3 extremely bad integrand behaviour occurs
!                             at some points of the integration
!                             interval.
!                         = 6 the input is invalid, because
!                             (epsabs.le.0 and
!                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
!                             result, abserr, neval, last, rlist(1) ,
!                             elist(1) and iord(1) are set to zero.
!                             alist(1) and blist(1) are set to a and b
!                             respectively.
!
!            alist   - double precision
!                      vector of dimension at least limit, the first
!                       last  elements of which are the left
!                      end points of the subintervals in the partition
!                      of the given integration range (a,b)
!
!            blist   - double precision
!                      vector of dimension at least limit, the first
!                       last  elements of which are the right
!                      end points of the subintervals in the partition
!                      of the given integration range (a,b)
!
!            rlist   - double precision
!                      vector of dimension at least limit, the first
!                       last  elements of which are the
!                      integral approximations on the subintervals
!
!            elist   - double precision
!                      vector of dimension at least limit, the first
!                       last  elements of which are the moduli of the
!                      absolute error estimates on the subintervals
!
!            iord    - integer
!                      vector of dimension at least limit, the first k
!                      elements of which are pointers to the
!                      error estimates over the subintervals,
!                      such that elist(iord(1)), ...,
!                      elist(iord(k)) form a decreasing sequence,
!                      with k = last if last.le.(limit/2+2), and
!                      k = limit+1-last otherwise
!
!            last    - integer
!                      number of subintervals actually produced in the
!                      subdivision process
!
!***references  (none)
!***routines called  d1mach,dqk15,dqk21,dqk31,
!                    dqk41,dqk51,dqk61,dqpsrt
!***end prologue  dqage
!
      double precision a,abserr,alist,area,area1,area12,area2,a1,a2,b, &
       blist,b1,b2,dabs,defabs,defab1,defab2,dmax1,elist,epmach, &
       epsabs,epsrel,errbnd,errmax,error1,error2,erro12,errsum, &
       resabs,result,rlist,uflow
      integer ier,iord,iroff1,iroff2,k,key,keyf,last,limit,maxerr,neval, &
       nrmax
!
      dimension alist(limit),blist(limit),elist(limit),iord(limit), &
       rlist(limit)
!
      double precision,external :: f,g,h
      

!
!            list of major variables
!            -----------------------
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                      (alist(i),blist(i))
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest
!                       error estimate
!           errmax    - elist(maxerr)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left subinterval
!           *****2    - variable for the right subinterval
!           last      - index for subdivision
!
!
!           machine dependent constants
!           ---------------------------
!
!           epmach  is the largest relative spacing.
!           uflow  is the smallest positive magnitude.
!
!***first executable statement  dqage
      epmach = d1mach(4)
      uflow = d1mach(1)
!
!           test on validity of parameters
!           ------------------------------
!
      ier = 0
      neval = 0
      last = 0
      result = 0.0d+00
      abserr = 0.0d+00
      alist(1) = a
      blist(1) = b
      rlist(1) = 0.0d+00
      elist(1) = 0.0d+00
      iord(1) = 0
      if(epsabs.le.0.0d+00.and. &
       epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) ier = 6 
      if(ier.eq.6) go to 999
!
!           first approximation to the integral
!           -----------------------------------
!
      keyf = key
      if(key.le.0) keyf = 1
      if(key.ge.7) keyf = 6
      neval = 0
      if(keyf.eq.1) then
        call dqk15_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      if(keyf.eq.2) then
        call dqk21_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      if(keyf.eq.3) then
        call dqk31_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      if(keyf.eq.4) then
        call dqk41_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      if(keyf.eq.5) then
        call dqk51_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      if(keyf.eq.6) then
        call dqk61_2d_outer(f,a,b,g,h,result,abserr,defabs,resabs, &
        epsabs,epsrel,key,limit)
      end if
      last = 1
      rlist(1) = result
      elist(1) = abserr
      iord(1) = 1
!
!           test on accuracy.
!
      errbnd = dmax1(epsabs,epsrel*dabs(result))
      if(abserr.le.0.5d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2
      if(limit.eq.1) ier = 1
      if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs) &
       .or.abserr.eq.0.0d+00) go to 60
!
!           initialization
!           --------------
!
!
      errmax = abserr
      maxerr = 1
      area = result
      errsum = abserr
      nrmax = 1
      iroff1 = 0
      iroff2 = 0
!
!           main do-loop
!           ------------
!
      do 30 last = 2,limit
!
!           bisect the subinterval with the largest error estimate.
!
        a1 = alist(maxerr)
        b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
        a2 = b1
        b2 = blist(maxerr)
        if(keyf.eq.1) then
            call dqk15_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.2) then
            call dqk21_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.3) then
            call dqk31_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.4) then
            call dqk41_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.5) then
            call dqk51_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.6) then
            call dqk61_2d_outer(f,a1,b1,g,h,area1,error1,resabs,defab1, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.1) then
            call dqk15_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.2) then
            call dqk21_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.3) then
            call dqk31_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.4) then
            call dqk41_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.5) then
            call dqk51_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
        if(keyf.eq.6) then
            call dqk61_2d_outer(f,a2,b2,g,h,area2,error2,resabs,defab2, &
        epsabs,epsrel,key,limit)
        end if
!
!           improve previous approximations to integral
!           and error and test for accuracy.
!
        neval = neval+1
        area12 = area1+area2
        erro12 = error1+error2
        errsum = errsum+erro12-errmax
        area = area+area12-rlist(maxerr)
        if(defab1.eq.error1.or.defab2.eq.error2) go to 5
        if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12) &
       .and.erro12.ge.0.99d+00*errmax) iroff1 = iroff1+1
        if(last.gt.10.and.erro12.gt.errmax) iroff2 = iroff2+1
    5   rlist(maxerr) = area1
        rlist(last) = area2
        errbnd = dmax1(epsabs,epsrel*dabs(area))
        if(errsum.le.errbnd) go to 8
!
!           test for roundoff error and eventually set error flag.
!
        if(iroff1.ge.6.or.iroff2.ge.20) ier = 2
!
!           set error flag in the case that the number of subintervals
!           equals limit.
!
        if(last.eq.limit) ier = 1
!
!           set error flag in the case of bad integrand behaviour
!           at a point of the integration range.
!
        if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03* &
       epmach)*(dabs(a2)+0.1d+04*uflow)) ier = 3
!
!           append the newly-created intervals to the list.
!
    8   if(error2.gt.error1) go to 10
        alist(last) = a2
        blist(maxerr) = b1
        blist(last) = b2
        elist(maxerr) = error1
        elist(last) = error2
        go to 20
   10   alist(maxerr) = a2
        alist(last) = a1
        blist(last) = b1
        rlist(maxerr) = area2
        rlist(last) = area1
        elist(maxerr) = error2
        elist(last) = error1
!
!           call subroutine dqpsrt to maintain the descending ordering
!           in the list of error estimates and select the subinterval
!           with the largest error estimate (to be bisected next).
!
   20   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
! ***jump out of do-loop
        if(ier.ne.0.or.errsum.le.errbnd) go to 40
   30 continue
!
!           compute final result.
!           ---------------------
!
   40 result = 0.0d+00
      do 50 k=1,last
        result = result+rlist(k)
   50 continue
      abserr = errsum
   60 if(keyf.ne.1) neval = (10*keyf+1)*(2*neval+1)
      if(keyf.eq.1) neval = 30*neval+15
  999 return
      end subroutine