!   fvn comment :
!   Modified version of the dqk61 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 dqk61_2d_inner(f,x,a,b,result,abserr,resabs,resasc)
!***begin prologue  dqk61
!***date written   800101   (yymmdd)
!***revision date  830518   (yymmdd)
!***category no.  h2a1a2
!***keywords  61-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 dabs(f) over (a,b)
!***description
!
!        integration rule
!        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 calling 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 61-point
!                    kronrod rule (resk) obtained by optimal addition of
!                    abscissae to the 30-point gauss rule (resg).
!
!           abserr - double precision
!                    estimate of the modulus of the absolute error,
!                    which should equal or exceed dabs(i-result)
!
!           resabs - double precision
!                    approximation to the integral j
!
!           resasc - double precision
!                    approximation to the integral of dabs(f-i/(b-a))
!
!
!***references  (none)
!***routines called  d1mach
!***end prologue  dqk61
!
      double precision a,dabsc,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(30),fv2(30),xgk(31),wgk(31),wg(15)
!
!           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 61-point kronrod rule
!                   xgk(2), xgk(4)  ... abscissae of the 30-point
!                   gauss rule
!                   xgk(1), xgk(3)  ... optimally added abscissae
!                   to the 30-point gauss rule
!
!           wgk   - weights of the 61-point kronrod rule
!
!           wg    - weigths of the 30-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.007968192496166605615465883474674d0 /
      data wg  (  2) / 0.018466468311090959142302131912047d0 /
      data wg  (  3) / 0.028784707883323369349719179611292d0 /
      data wg  (  4) / 0.038799192569627049596801936446348d0 /
      data wg  (  5) / 0.048402672830594052902938140422808d0 /
      data wg  (  6) / 0.057493156217619066481721689402056d0 /
      data wg  (  7) / 0.065974229882180495128128515115962d0 /
      data wg  (  8) / 0.073755974737705206268243850022191d0 /
      data wg  (  9) / 0.080755895229420215354694938460530d0 /
      data wg  ( 10) / 0.086899787201082979802387530715126d0 /
      data wg  ( 11) / 0.092122522237786128717632707087619d0 /
      data wg  ( 12) / 0.096368737174644259639468626351810d0 /
      data wg  ( 13) / 0.099593420586795267062780282103569d0 /
      data wg  ( 14) / 0.101762389748405504596428952168554d0 /
      data wg  ( 15) / 0.102852652893558840341285636705415d0 /
!
      data xgk (  1) / 0.999484410050490637571325895705811d0 /
      data xgk (  2) / 0.996893484074649540271630050918695d0 /
      data xgk (  3) / 0.991630996870404594858628366109486d0 /
      data xgk (  4) / 0.983668123279747209970032581605663d0 /
      data xgk (  5) / 0.973116322501126268374693868423707d0 /
      data xgk (  6) / 0.960021864968307512216871025581798d0 /
      data xgk (  7) / 0.944374444748559979415831324037439d0 /
      data xgk (  8) / 0.926200047429274325879324277080474d0 /
      data xgk (  9) / 0.905573307699907798546522558925958d0 /
      data xgk ( 10) / 0.882560535792052681543116462530226d0 /
      data xgk ( 11) / 0.857205233546061098958658510658944d0 /
      data xgk ( 12) / 0.829565762382768397442898119732502d0 /
      data xgk ( 13) / 0.799727835821839083013668942322683d0 /
      data xgk ( 14) / 0.767777432104826194917977340974503d0 /
      data xgk ( 15) / 0.733790062453226804726171131369528d0 /
      data xgk ( 16) / 0.697850494793315796932292388026640d0 /
      data xgk ( 17) / 0.660061064126626961370053668149271d0 /
      data xgk ( 18) / 0.620526182989242861140477556431189d0 /
      data xgk ( 19) / 0.579345235826361691756024932172540d0 /
      data xgk ( 20) / 0.536624148142019899264169793311073d0 /
      data xgk ( 21) / 0.492480467861778574993693061207709d0 /
      data xgk ( 22) / 0.447033769538089176780609900322854d0 /
      data xgk ( 23) / 0.400401254830394392535476211542661d0 /
      data xgk ( 24) / 0.352704725530878113471037207089374d0 /
      data xgk ( 25) / 0.304073202273625077372677107199257d0 /
      data xgk ( 26) / 0.254636926167889846439805129817805d0 /
      data xgk ( 27) / 0.204525116682309891438957671002025d0 /
      data xgk ( 28) / 0.153869913608583546963794672743256d0 /
      data xgk ( 29) / 0.102806937966737030147096751318001d0 /
      data xgk ( 30) / 0.051471842555317695833025213166723d0 /
      data xgk ( 31) / 0.000000000000000000000000000000000d0 /
!
      data wgk (  1) / 0.001389013698677007624551591226760d0 /
      data wgk (  2) / 0.003890461127099884051267201844516d0 /
      data wgk (  3) / 0.006630703915931292173319826369750d0 /
      data wgk (  4) / 0.009273279659517763428441146892024d0 /
      data wgk (  5) / 0.011823015253496341742232898853251d0 /
      data wgk (  6) / 0.014369729507045804812451432443580d0 /
      data wgk (  7) / 0.016920889189053272627572289420322d0 /
      data wgk (  8) / 0.019414141193942381173408951050128d0 /
      data wgk (  9) / 0.021828035821609192297167485738339d0 /
      data wgk ( 10) / 0.024191162078080601365686370725232d0 /
      data wgk ( 11) / 0.026509954882333101610601709335075d0 /
      data wgk ( 12) / 0.028754048765041292843978785354334d0 /
      data wgk ( 13) / 0.030907257562387762472884252943092d0 /
      data wgk ( 14) / 0.032981447057483726031814191016854d0 /
      data wgk ( 15) / 0.034979338028060024137499670731468d0 /
      data wgk ( 16) / 0.036882364651821229223911065617136d0 /
      data wgk ( 17) / 0.038678945624727592950348651532281d0 /
      data wgk ( 18) / 0.040374538951535959111995279752468d0 /
      data wgk ( 19) / 0.041969810215164246147147541285970d0 /
      data wgk ( 20) / 0.043452539701356069316831728117073d0 /
      data wgk ( 21) / 0.044814800133162663192355551616723d0 /
      data wgk ( 22) / 0.046059238271006988116271735559374d0 /
      data wgk ( 23) / 0.047185546569299153945261478181099d0 /
      data wgk ( 24) / 0.048185861757087129140779492298305d0 /
      data wgk ( 25) / 0.049055434555029778887528165367238d0 /
      data wgk ( 26) / 0.049795683427074206357811569379942d0 /
      data wgk ( 27) / 0.050405921402782346840893085653585d0 /
      data wgk ( 28) / 0.050881795898749606492297473049805d0 /
      data wgk ( 29) / 0.051221547849258772170656282604944d0 /
      data wgk ( 30) / 0.051426128537459025933862879215781d0 /
      data wgk ( 31) / 0.051494729429451567558340433647099d0 /
!
!           list of major variables
!           -----------------------
!
!           centr  - mid point of the interval
!           hlgth  - half-length of the interval
!           dabsc  - abscissa
!           fval*  - function value
!           resg   - result of the 30-point gauss rule
!           resk   - result of the 61-point kronrod rule
!           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.
!
      epmach = d1mach(4)
      uflow = d1mach(1)
!
      centr = 0.5d+00*(b+a)
      hlgth = 0.5d+00*(b-a)
      dhlgth = dabs(hlgth)
!
!           compute the 61-point kronrod approximation to the
!           integral, and estimate the absolute error.
!
!***first executable statement  dqk61
      resg = 0.0d+00
      fc = f(x,centr)
      resk = wgk(31)*fc
      resabs = dabs(resk)
      do 10 j=1,15
        jtw = j*2
        dabsc = hlgth*xgk(jtw)
        fval1 = f(x,centr-dabsc)
        fval2 = f(x,centr+dabsc)
        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,15
        jtwm1 = j*2-1
        dabsc = hlgth*xgk(jtwm1)
        fval1 = f(x,centr-dabsc)
        fval2 = f(x,centr+dabsc)
        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(31)*dabs(fc-reskh)
      do 20 j=1,30
        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