!   fvn comment :
!   Modified version of the dqk21 quadpack routine from http://www.netlib.org/quadpack
!
!   + The external 'f' function is a 2 parameters function f(x,y). The routine
!   takes one more parameter 'x' and evaluate the integral of f against y between a and b
!   for a given x
      subroutine dqk21_2d_inner(f,x,a,b,result,abserr,resabs,resasc)
!***begin prologue  dqk21
!***date written   800101   (yymmdd)
!***revision date  830518   (yymmdd)
!***category no.  h2a1a2
!***keywords  21-point gauss-kronrod rules
!***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
!           de doncker,elise,appl. math. & progr. div. - k.u.leuven
!***purpose  to compute i = integral of f over (a,b), with error
!                           estimate
!                       j = integral of abs(f) over (a,b)
!***description
!
!           integration rules
!           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
!
!            on return
!              result - double precision
!                       approximation to the integral i
!                       result is computed by applying the 21-point
!                       kronrod rule (resk) obtained by optimal addition
!                       of abscissae to the 10-point gauss rule (resg).
!
!              abserr - double precision
!                       estimate of the modulus of the absolute error,
!                       which should not exceed abs(i-result)
!
!              resabs - double precision
!                       approximation to the integral j
!
!              resasc - double precision
!                       approximation to the integral of abs(f-i/(b-a))
!                       over (a,b)
!
!***references  (none)
!***routines called  d1mach
!***end prologue  dqk21
!
      double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1, &
       epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc, &
       resg,resk,reskh,result,uflow,wg,wgk,xgk,x
      integer j,jtw,jtwm1
      external f
!
      dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11)
!
!           the abscissae and weights are given for the interval (-1,1).
!           because of symmetry only the positive abscissae and their
!           corresponding weights are given.
!
!           xgk    - abscissae of the 21-point kronrod rule
!                    xgk(2), xgk(4), ...  abscissae of the 10-point
!                    gauss rule
!                    xgk(1), xgk(3), ...  abscissae which are optimally
!                    added to the 10-point gauss rule
!
!           wgk    - weights of the 21-point kronrod rule
!
!           wg     - weights of the 10-point gauss rule
!
!
! gauss quadrature weights and kronron quadrature abscissae and weights
! as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
! bell labs, nov. 1981.
!
      data wg  (  1) / 0.066671344308688137593568809893332d0 /
      data wg  (  2) / 0.149451349150580593145776339657697d0 /
      data wg  (  3) / 0.219086362515982043995534934228163d0 /
      data wg  (  4) / 0.269266719309996355091226921569469d0 /
      data wg  (  5) / 0.295524224714752870173892994651338d0 /
!
      data xgk (  1) / 0.995657163025808080735527280689003d0 /
      data xgk (  2) / 0.973906528517171720077964012084452d0 /
      data xgk (  3) / 0.930157491355708226001207180059508d0 /
      data xgk (  4) / 0.865063366688984510732096688423493d0 /
      data xgk (  5) / 0.780817726586416897063717578345042d0 /
      data xgk (  6) / 0.679409568299024406234327365114874d0 /
      data xgk (  7) / 0.562757134668604683339000099272694d0 /
      data xgk (  8) / 0.433395394129247190799265943165784d0 /
      data xgk (  9) / 0.294392862701460198131126603103866d0 /
      data xgk ( 10) / 0.148874338981631210884826001129720d0 /
      data xgk ( 11) / 0.000000000000000000000000000000000d0 /
!
      data wgk (  1) / 0.011694638867371874278064396062192d0 /
      data wgk (  2) / 0.032558162307964727478818972459390d0 /
      data wgk (  3) / 0.054755896574351996031381300244580d0 /
      data wgk (  4) / 0.075039674810919952767043140916190d0 /
      data wgk (  5) / 0.093125454583697605535065465083366d0 /
      data wgk (  6) / 0.109387158802297641899210590325805d0 /
      data wgk (  7) / 0.123491976262065851077958109831074d0 /
      data wgk (  8) / 0.134709217311473325928054001771707d0 /
      data wgk (  9) / 0.142775938577060080797094273138717d0 /
      data wgk ( 10) / 0.147739104901338491374841515972068d0 /
      data wgk ( 11) / 0.149445554002916905664936468389821d0 /
!
!
!           list of major variables
!           -----------------------
!
!           centr  - mid point of the interval
!           hlgth  - half-length of the interval
!           absc   - abscissa
!           fval*  - function value
!           resg   - result of the 10-point gauss formula
!           resk   - result of the 21-point kronrod formula
!           reskh  - approximation to the mean value of f over (a,b),
!                    i.e. to i/(b-a)
!
!
!           machine dependent constants
!           ---------------------------
!
!           epmach is the largest relative spacing.
!           uflow is the smallest positive magnitude.
!
!***first executable statement  dqk21
      epmach = d1mach(4)
      uflow = d1mach(1)
!
      centr = 0.5d+00*(a+b)
      hlgth = 0.5d+00*(b-a)
      dhlgth = dabs(hlgth)
!
!           compute the 21-point kronrod approximation to
!           the integral, and estimate the absolute error.
!
      resg = 0.0d+00
      fc = f(x,centr)
      resk = wgk(11)*fc
      resabs = dabs(resk)
      do 10 j=1,5
        jtw = 2*j
        absc = hlgth*xgk(jtw)
        fval1 = f(x,centr-absc)
        fval2 = f(x,centr+absc)
        fv1(jtw) = fval1
        fv2(jtw) = fval2
        fsum = fval1+fval2
        resg = resg+wg(j)*fsum
        resk = resk+wgk(jtw)*fsum
        resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
   10 continue
      do 15 j = 1,5
        jtwm1 = 2*j-1
        absc = hlgth*xgk(jtwm1)
        fval1 = f(x,centr-absc)
        fval2 = f(x,centr+absc)
        fv1(jtwm1) = fval1
        fv2(jtwm1) = fval2
        fsum = fval1+fval2
        resk = resk+wgk(jtwm1)*fsum
        resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
   15 continue
      reskh = resk*0.5d+00
      resasc = wgk(11)*dabs(fc-reskh)
      do 20 j=1,10
        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
   20 continue
      result = resk*hlgth
      resabs = resabs*dhlgth
      resasc = resasc*dhlgth
      abserr = dabs((resk-resg)*hlgth)
      if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) &
       abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
      if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1 &
       ((epmach*0.5d+02)*resabs,abserr)
      return
      end subroutine