print in fort.33 also for tau>taucr; changed test for Npl>1 in after_onestep

Makefile

@@ -11,7 +11,7 @@ vpath %.f src
 
 # Fortran compiler name and flags
 FC=gfortran
-FFLAGS=-O3 
+FFLAGS=-O3
 
 DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
 

src/gray.f

@@ -227,13 +227,13 @@ c
           call gwork(j,k)
 c 
           if(ierr.gt.0) then
-!            print*,'    IERR = ', ierr
+            print*,'    IERR = ', ierr
             if(ierr.eq.97) then
 c              igrad=0
 c              ierr=0 
             else  
               istop=1
-              exit
+c              exit
             end if
           end if
 
@@ -353,6 +353,7 @@ c        print*,'istep = ',i
 c
       if (istpl.eq.istpl0) istpl=0
       istpl=istpl+1
+      
       return
       end
 c
@@ -482,7 +483,8 @@ c     central ray only end
 
 c     print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
       if(istpl.eq.istpl0)  then
-        if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
+        if(j.eq.nrayr) then
+c        if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
           write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
      .          psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
         end if
@@ -4967,8 +4969,6 @@ c              call dqagi(fhermit,bound,2,epsa,epsr,resfh,
             end if
           end do
         end do
-c      write(83,'(12(1x,e12.5))') 
-c     .       yg,rr(1,0,1),rr(1,1,1),rr(1,2,1)
       
       if(iwarm.gt.10) return
 c

src/grayl.f

@@ -10902,3 +10902,780 @@ c
       errest = 2.0d0*errest
       go to 82
       end
+
+      subroutine dqagp(f,a,b,npts2,points,epsabs,epsrel,result,abserr,
+     *   neval,ier,leniw,lenw,last,iwork,work)
+c***begin prologue  dqagp
+c***date written   800101   (yymmdd)
+c***revision date  830518   (yymmdd)
+c***category no.  h2a2a1
+c***keywords  automatic integrator, general-purpose,
+c             singularities at user specified points,
+c             extrapolation, globally adaptive
+c***author  piessens,robert,appl. math. & progr. div - k.u.leuven
+c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
+c***purpose  the routine calculates an approximation result to a given
+c            definite integral i = integral of f over (a,b),
+c            hopefully satisfying following claim for accuracy
+c            break points of the integration interval, where local
+c            difficulties of the integrand may occur (e.g.
+c            singularities, discontinuities), are provided by the user.
+c***description
+c
+c        computation of a definite integral
+c        standard fortran subroutine
+c        double precision version
+c
+c        parameters
+c         on entry
+c            f      - double precision
+c                     function subprogram defining the integrand
+c                     function f(x). the actual name for f needs to be
+c                     declared e x t e r n a l in the driver program.
+c
+c            a      - double precision
+c                     lower limit of integration
+c
+c            b      - double precision
+c                     upper limit of integration
+c
+c            npts2  - integer
+c                     number equal to two more than the number of
+c                     user-supplied break points within the integration
+c                     range, npts.ge.2.
+c                     if npts2.lt.2, the routine will end with ier = 6.
+c
+c            points - double precision
+c                     vector of dimension npts2, the first (npts2-2)
+c                     elements of which are the user provided break
+c                     points. if these points do not constitute an
+c                     ascending sequence there will be an automatic
+c                     sorting.
+c
+c            epsabs - double precision
+c                     absolute accuracy requested
+c            epsrel - double precision
+c                     relative accuracy requested
+c                     if  epsabs.le.0
+c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+c                     the routine will end with ier = 6.
+c
+c         on return
+c            result - double precision
+c                     approximation to the integral
+c
+c            abserr - double precision
+c                     estimate of the modulus of the absolute error,
+c                     which should equal or exceed abs(i-result)
+c
+c            neval  - integer
+c                     number of integrand evaluations
+c
+c            ier    - integer
+c                     ier = 0 normal and reliable termination of the
+c                             routine. it is assumed that the requested
+c                             accuracy has been achieved.
+c                     ier.gt.0 abnormal termination of the routine.
+c                             the estimates for integral and error are
+c                             less reliable. it is assumed that the
+c                             requested accuracy has not been achieved.
+c            error messages
+c                     ier = 1 maximum number of subdivisions allowed
+c                             has been achieved. one can allow more
+c                             subdivisions by increasing the value of
+c                             limit (and taking the according dimension
+c                             adjustments into account). however, if
+c                             this yields no improvement it is advised
+c                             to analyze the integrand in order to
+c                             determine the integration difficulties. if
+c                             the position of a local difficulty can be
+c                             determined (i.e. singularity,
+c                             discontinuity within the interval), it
+c                             should be supplied to the routine as an
+c                             element of the vector points. if necessary
+c                             an appropriate special-purpose integrator
+c                             must be used, which is designed for
+c                             handling the type of difficulty involved.
+c                         = 2 the occurrence of roundoff error is
+c                             detected, which prevents the requested
+c                             tolerance from being achieved.
+c                             the error may be under-estimated.
+c                         = 3 extremely bad integrand behaviour occurs
+c                             at some points of the integration
+c                             interval.
+c                         = 4 the algorithm does not converge.
+c                             roundoff error is detected in the
+c                             extrapolation table.
+c                             it is presumed that the requested
+c                             tolerance cannot be achieved, and that
+c                             the returned result is the best which
+c                             can be obtained.
+c                         = 5 the integral is probably divergent, or
+c                             slowly convergent. it must be noted that
+c                             divergence can occur with any other value
+c                             of ier.gt.0.
+c                         = 6 the input is invalid because
+c                             npts2.lt.2 or
+c                             break points are specified outside
+c                             the integration range or
+c                             (epsabs.le.0 and
+c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
+c                             result, abserr, neval, last are set to
+c                             zero. exept when leniw or lenw or npts2 is
+c                             invalid, iwork(1), iwork(limit+1),
+c                             work(limit*2+1) and work(limit*3+1)
+c                             are set to zero.
+c                             work(1) is set to a and work(limit+1)
+c                             to b (where limit = (leniw-npts2)/2).
+c
+c         dimensioning parameters
+c            leniw - integer
+c                    dimensioning parameter for iwork
+c                    leniw determines limit = (leniw-npts2)/2,
+c                    which is the maximum number of subintervals in the
+c                    partition of the given integration interval (a,b),
+c                    leniw.ge.(3*npts2-2).
+c                    if leniw.lt.(3*npts2-2), the routine will end with
+c                    ier = 6.
+c
+c            lenw  - integer
+c                    dimensioning parameter for work
+c                    lenw must be at least leniw*2-npts2.
+c                    if lenw.lt.leniw*2-npts2, the routine will end
+c                    with ier = 6.
+c
+c            last  - integer
+c                    on return, last equals the number of subintervals
+c                    produced in the subdivision process, which
+c                    determines the number of significant elements
+c                    actually in the work arrays.
+c
+c         work arrays
+c            iwork - integer
+c                    vector of dimension at least leniw. on return,
+c                    the first k elements of which contain
+c                    pointers to the error estimates over the
+c                    subintervals, such that work(limit*3+iwork(1)),...,
+c                    work(limit*3+iwork(k)) form a decreasing
+c                    sequence, with k = last if last.le.(limit/2+2), and
+c                    k = limit+1-last otherwise
+c                    iwork(limit+1), ...,iwork(limit+last) contain the
+c                     subdivision levels of the subintervals, i.e.
+c                     if (aa,bb) is a subinterval of (p1,p2)
+c                     where p1 as well as p2 is a user-provided
+c                     break point or integration limit, then (aa,bb) has
+c                     level l if abs(bb-aa) = abs(p2-p1)*2**(-l),
+c                    iwork(limit*2+1), ..., iwork(limit*2+npts2) have
+c                     no significance for the user,
+c                    note that limit = (leniw-npts2)/2.
+c
+c            work  - double precision
+c                    vector of dimension at least lenw
+c                    on return
+c                    work(1), ..., work(last) contain the left
+c                     end points of the subintervals in the
+c                     partition of (a,b),
+c                    work(limit+1), ..., work(limit+last) contain
+c                     the right end points,
+c                    work(limit*2+1), ..., work(limit*2+last) contain
+c                     the integral approximations over the subintervals,
+c                    work(limit*3+1), ..., work(limit*3+last)
+c                     contain the corresponding error estimates,
+c                    work(limit*4+1), ..., work(limit*4+npts2)
+c                     contain the integration limits and the
+c                     break points sorted in an ascending sequence.
+c                    note that limit = (leniw-npts2)/2.
+c
+c***references  (none)
+c***routines called  dqagpe,xerror
+c***end prologue  dqagp
+c
+      double precision a,abserr,b,epsabs,epsrel,f,points,result,work
+      integer ier,iwork,last,leniw,lenw,limit,lvl,l1,l2,l3,l4,neval,
+     *  npts2
+c
+      dimension iwork(leniw),points(npts2),work(lenw)
+c
+      external f
+c
+c         check validity of limit and lenw.
+c
+c***first executable statement  dqagp
+      ier = 6
+      neval = 0
+      last = 0
+      result = 0.0d+00
+      abserr = 0.0d+00
+      if(leniw.lt.(3*npts2-2).or.lenw.lt.(leniw*2-npts2).or.npts2.lt.2)
+     *  go to 10
+c
+c         prepare call for dqagpe.
+c
+      limit = (leniw-npts2)/2
+      l1 = limit+1
+      l2 = limit+l1
+      l3 = limit+l2
+      l4 = limit+l3
+c
+      call dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,abserr,
+     *  neval,ier,work(1),work(l1),work(l2),work(l3),work(l4),
+     *  iwork(1),iwork(l1),iwork(l2),last)
+c
+c         call error handler if necessary.
+c
+      lvl = 0
+10    if(ier.eq.6) lvl = 1
+      if(ier.ne.0) print*,'habnormal return from dqaps',ier,lvl
+      return
+      end
+
+      subroutine dqagpe(f,a,b,npts2,points,epsabs,epsrel,limit,result,
+     *   abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,
+     *   last)
+c***begin prologue  dqagpe
+c***date written   800101   (yymmdd)
+c***revision date  830518   (yymmdd)
+c***category no.  h2a2a1
+c***keywords  automatic integrator, general-purpose,
+c             singularities at user specified points,
+c             extrapolation, globally adaptive.
+c***author  piessens,robert ,appl. math. & progr. div. - k.u.leuven
+c           de doncker,elise,appl. math. & progr. div. - k.u.leuven
+c***purpose  the routine calculates an approximation result to a given
+c            definite integral i = integral of f over (a,b), hopefully
+c            satisfying following claim for accuracy abs(i-result).le.
+c            max(epsabs,epsrel*abs(i)). break points of the integration
+c            interval, where local difficulties of the integrand may
+c            occur(e.g. singularities,discontinuities),provided by user.
+c***description
+c
+c        computation of a definite integral
+c        standard fortran subroutine
+c        double precision version
+c
+c        parameters
+c         on entry
+c            f      - double precision
+c                     function subprogram defining the integrand
+c                     function f(x). the actual name for f needs to be
+c                     declared e x t e r n a l in the driver program.
+c
+c            a      - double precision
+c                     lower limit of integration
+c
+c            b      - double precision
+c                     upper limit of integration
+c
+c            npts2  - integer
+c                     number equal to two more than the number of
+c                     user-supplied break points within the integration
+c                     range, npts2.ge.2.
+c                     if npts2.lt.2, the routine will end with ier = 6.
+c
+c            points - double precision
+c                     vector of dimension npts2, the first (npts2-2)
+c                     elements of which are the user provided break
+c                     points. if these points do not constitute an
+c                     ascending sequence there will be an automatic
+c                     sorting.
+c
+c            epsabs - double precision
+c                     absolute accuracy requested
+c            epsrel - double precision
+c                     relative accuracy requested
+c                     if  epsabs.le.0
+c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
+c                     the routine will end with ier = 6.
+c
+c            limit  - integer
+c                     gives an upper bound on the number of subintervals
+c                     in the partition of (a,b), limit.ge.npts2
+c                     if limit.lt.npts2, the routine will end with
+c                     ier = 6.
+c
+c         on return
+c            result - double precision
+c                     approximation to the integral
+c
+c            abserr - double precision
+c                     estimate of the modulus of the absolute error,
+c                     which should equal or exceed abs(i-result)
+c
+c            neval  - integer
+c                     number of integrand evaluations
+c
+c            ier    - integer
+c                     ier = 0 normal and reliable termination of the
+c                             routine. it is assumed that the requested
+c                             accuracy has been achieved.
+c                     ier.gt.0 abnormal termination of the routine.
+c                             the estimates for integral and error are
+c                             less reliable. it is assumed that the
+c                             requested accuracy has not been achieved.
+c            error messages
+c                     ier = 1 maximum number of subdivisions allowed
+c                             has been achieved. one can allow more
+c                             subdivisions by increasing the value of
+c                             limit (and taking the according dimension
+c                             adjustments into account). however, if
+c                             this yields no improvement it is advised
+c                             to analyze the integrand in order to
+c                             determine the integration difficulties. if
+c                             the position of a local difficulty can be
+c                             determined (i.e. singularity,
+c                             discontinuity within the interval), it
+c                             should be supplied to the routine as an
+c                             element of the vector points. if necessary
+c                             an appropriate special-purpose integrator
+c                             must be used, which is designed for
+c                             handling the type of difficulty involved.
+c                         = 2 the occurrence of roundoff error is
+c                             detected, which prevents the requested
+c                             tolerance from being achieved.
+c                             the error may be under-estimated.
+c                         = 3 extremely bad integrand behaviour occurs
+c                             at some points of the integration
+c                             interval.
+c                         = 4 the algorithm does not converge.
+c                             roundoff error is detected in the
+c                             extrapolation table. it is presumed that
+c                             the requested tolerance cannot be
+c                             achieved, and that the returned result is
+c                             the best which can be obtained.
+c                         = 5 the integral is probably divergent, or
+c                             slowly convergent. it must be noted that
+c                             divergence can occur with any other value
+c                             of ier.gt.0.
+c                         = 6 the input is invalid because
+c                             npts2.lt.2 or
+c                             break points are specified outside
+c                             the integration range or
+c                             (epsabs.le.0 and
+c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
+c                             or limit.lt.npts2.
+c                             result, abserr, neval, last, rlist(1),
+c                             and elist(1) are set to zero. alist(1) and
+c                             blist(1) are set to a and b respectively.
+c
+c            alist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the left end points
+c                     of the subintervals in the partition of the given
+c                     integration range (a,b)
+c
+c            blist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the right end points
+c                     of the subintervals in the partition of the given
+c                     integration range (a,b)
+c
+c            rlist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the integral
+c                     approximations on the subintervals
+c
+c            elist  - double precision
+c                     vector of dimension at least limit, the first
+c                      last  elements of which are the moduli of the
+c                     absolute error estimates on the subintervals
+c
+c            pts    - double precision
+c                     vector of dimension at least npts2, containing the
+c                     integration limits and the break points of the
+c                     interval in ascending sequence.
+c
+c            level  - integer
+c                     vector of dimension at least limit, containing the
+c                     subdivision levels of the subinterval, i.e. if
+c                     (aa,bb) is a subinterval of (p1,p2) where p1 as
+c                     well as p2 is a user-provided break point or
+c                     integration limit, then (aa,bb) has level l if
+c                     abs(bb-aa) = abs(p2-p1)*2**(-l).
+c
+c            ndin   - integer
+c                     vector of dimension at least npts2, after first
+c                     integration over the intervals (pts(i)),pts(i+1),
+c                     i = 0,1, ..., npts2-2, the error estimates over
+c                     some of the intervals may have been increased
+c                     artificially, in order to put their subdivision
+c                     forward. if this happens for the subinterval
+c                     numbered k, ndin(k) is put to 1, otherwise
+c                     ndin(k) = 0.
+c
+c            iord   - integer
+c                     vector of dimension at least limit, the first k
+c                     elements of which are pointers to the
+c                     error estimates over the subintervals,
+c                     such that elist(iord(1)), ..., elist(iord(k))
+c                     form a decreasing sequence, with k = last
+c                     if last.le.(limit/2+2), and k = limit+1-last
+c                     otherwise
+c
+c            last   - integer
+c                     number of subintervals actually produced in the
+c                     subdivisions process
+c
+c***references  (none)
+c***routines called  d1mach,dqelg,dqk21,dqpsrt
+c***end prologue  dqagpe
+      double precision a,abseps,abserr,alist,area,area1,area12,area2,a1,
+     *  a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,dmax1,dmin1,
+     *  dres,d1mach,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,
+     *  errmax,error1,erro12,error2,errsum,ertest,f,oflow,points,pts,
+     *  resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp,uflow
+      integer i,id,ier,ierro,ind1,ind2,iord,ip1,iroff1,iroff2,iroff3,j,
+     *  jlow,jupbnd,k,ksgn,ktmin,last,levcur,level,levmax,limit,maxerr,
+     *  ndin,neval,nint,nintp1,npts,npts2,nres,nrmax,numrl2
+      logical extrap,noext
+c
+c
+      dimension alist(limit),blist(limit),elist(limit),iord(limit),
+     *  level(limit),ndin(npts2),points(npts2),pts(npts2),res3la(3),
+     *  rlist(limit),rlist2(52)
+c
+      external f
+c
+c            the dimension of rlist2 is determined by the value of
+c            limexp in subroutine epsalg (rlist2 should be of dimension
+c            (limexp+2) at least).
+c
+c
+c            list of major variables
+c            -----------------------
+c
+c           alist     - list of left end points of all subintervals
+c                       considered up to now
+c           blist     - list of right end points of all subintervals
+c                       considered up to now
+c           rlist(i)  - approximation to the integral over
+c                       (alist(i),blist(i))
+c           rlist2    - array of dimension at least limexp+2
+c                       containing the part of the epsilon table which
+c                       is still needed for further computations
+c           elist(i)  - error estimate applying to rlist(i)
+c           maxerr    - pointer to the interval with largest error
+c                       estimate
+c           errmax    - elist(maxerr)
+c           erlast    - error on the interval currently subdivided
+c                       (before that subdivision has taken place)
+c           area      - sum of the integrals over the subintervals
+c           errsum    - sum of the errors over the subintervals
+c           errbnd    - requested accuracy max(epsabs,epsrel*
+c                       abs(result))
+c           *****1    - variable for the left subinterval
+c           *****2    - variable for the right subinterval
+c           last      - index for subdivision
+c           nres      - number of calls to the extrapolation routine
+c           numrl2    - number of elements in rlist2. if an appropriate
+c                       approximation to the compounded integral has
+c                       been obtained, it is put in rlist2(numrl2) after
+c                       numrl2 has been increased by one.
+c           erlarg    - sum of the errors over the intervals larger
+c                       than the smallest interval considered up to now
+c           extrap    - logical variable denoting that the routine
+c                       is attempting to perform extrapolation. i.e.
+c                       before subdividing the smallest interval we
+c                       try to decrease the value of erlarg.
+c           noext     - logical variable denoting that extrapolation is
+c                       no longer allowed (true-value)
+c
+c            machine dependent constants
+c            ---------------------------
+c
+c           epmach is the largest relative spacing.
+c           uflow is the smallest positive magnitude.
+c           oflow is the largest positive magnitude.
+c
+c***first executable statement  dqagpe
+      epmach = d1mach(4)
+c
+c            test on validity of parameters
+c            -----------------------------
+c
+      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
+      level(1) = 0
+      npts = npts2-2
+      if(npts2.lt.2.or.limit.le.npts.or.(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
+c
+c            if any break points are provided, sort them into an
+c            ascending sequence.
+c
+      sign = 1.0d+00
+      if(a.gt.b) sign = -1.0d+00
+      pts(1) = dmin1(a,b)
+      if(npts.eq.0) go to 15
+      do 10 i = 1,npts
+        pts(i+1) = points(i)
+   10 continue
+   15 pts(npts+2) = dmax1(a,b)
+      nint = npts+1
+      a1 = pts(1)
+      if(npts.eq.0) go to 40
+      nintp1 = nint+1
+      do 20 i = 1,nint
+        ip1 = i+1
+        do 20 j = ip1,nintp1
+          if(pts(i).le.pts(j)) go to 20
+          temp = pts(i)
+          pts(i) = pts(j)
+          pts(j) = temp
+   20 continue
+      if(pts(1).ne.dmin1(a,b).or.pts(nintp1).ne.dmax1(a,b)) ier = 6
+      if(ier.eq.6) go to 999
+c
+c            compute first integral and error approximations.
+c            ------------------------------------------------
+c
+   40 resabs = 0.0d+00
+      do 50 i = 1,nint
+        b1 = pts(i+1)
+        call dqk21(f,a1,b1,area1,error1,defabs,resa)
+        abserr = abserr+error1
+        result = result+area1
+        ndin(i) = 0
+        if(error1.eq.resa.and.error1.ne.0.0d+00) ndin(i) = 1
+        resabs = resabs+defabs
+        level(i) = 0
+        elist(i) = error1
+        alist(i) = a1
+        blist(i) = b1
+        rlist(i) = area1
+        iord(i) = i
+        a1 = b1
+   50 continue
+      errsum = 0.0d+00
+      do 55 i = 1,nint
+        if(ndin(i).eq.1) elist(i) = abserr
+        errsum = errsum+elist(i)
+   55 continue
+c
+c           test on accuracy.
+c
+      last = nint
+      neval = 21*nint
+      dres = dabs(result)
+      errbnd = dmax1(epsabs,epsrel*dres)
+      if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2
+      if(nint.eq.1) go to 80
+      do 70 i = 1,npts
+        jlow = i+1
+        ind1 = iord(i)
+        do 60 j = jlow,nint
+          ind2 = iord(j)
+          if(elist(ind1).gt.elist(ind2)) go to 60
+          ind1 = ind2
+          k = j
+   60   continue
+        if(ind1.eq.iord(i)) go to 70
+        iord(k) = iord(i)
+        iord(i) = ind1
+   70 continue
+      if(limit.lt.npts2) ier = 1
+   80 if(ier.ne.0.or.abserr.le.errbnd) go to 210
+c
+c           initialization
+c           --------------
+c
+      rlist2(1) = result
+      maxerr = iord(1)
+      errmax = elist(maxerr)
+      area = result
+      nrmax = 1
+      nres = 0
+      numrl2 = 1
+      ktmin = 0
+      extrap = .false.
+      noext = .false.
+      erlarg = errsum
+      ertest = errbnd
+      levmax = 1
+      iroff1 = 0
+      iroff2 = 0
+      iroff3 = 0
+      ierro = 0
+      uflow = d1mach(1)
+      oflow = d1mach(2)
+      abserr = oflow
+      ksgn = -1
+      if(dres.ge.(0.1d+01-0.5d+02*epmach)*resabs) ksgn = 1
+c
+c           main do-loop
+c           ------------
+c
+      do 160 last = npts2,limit
+c
+c           bisect the subinterval with the nrmax-th largest error
+c           estimate.
+c
+        levcur = level(maxerr)+1
+        a1 = alist(maxerr)
+        b1 = 0.5d+00*(alist(maxerr)+blist(maxerr))
+        a2 = b1
+        b2 = blist(maxerr)
+        erlast = errmax
+        call dqk21(f,a1,b1,area1,error1,resa,defab1)
+        call dqk21(f,a2,b2,area2,error2,resa,defab2)
+c
+c           improve previous approximations to integral
+c           and error and test for accuracy.
+c
+        neval = neval+42
+        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 95
+        if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12)
+     *  .or.erro12.lt.0.99d+00*errmax) go to 90
+        if(extrap) iroff2 = iroff2+1
+        if(.not.extrap) iroff1 = iroff1+1
+   90   if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1
+   95   level(maxerr) = levcur
+        level(last) = levcur
+        rlist(maxerr) = area1
+        rlist(last) = area2
+        errbnd = dmax1(epsabs,epsrel*dabs(area))
+c
+c           test for roundoff error and eventually set error flag.
+c
+        if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2
+        if(iroff2.ge.5) ierro = 3
+c
+c           set error flag in the case that the number of
+c           subintervals equals limit.
+c
+        if(last.eq.limit) ier = 1
+c
+c           set error flag in the case of bad integrand behaviour
+c           at a point of the integration range
+c
+        if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)*
+     *  (dabs(a2)+0.1d+04*uflow)) ier = 4
+c
+c           append the newly-created intervals to the list.
+c
+        if(error2.gt.error1) go to 100
+        alist(last) = a2
+        blist(maxerr) = b1
+        blist(last) = b2
+        elist(maxerr) = error1
+        elist(last) = error2
+        go to 110
+  100   alist(maxerr) = a2
+        alist(last) = a1
+        blist(last) = b1
+        rlist(maxerr) = area2
+        rlist(last) = area1
+        elist(maxerr) = error2
+        elist(last) = error1
+c
+c           call subroutine dqpsrt to maintain the descending ordering
+c           in the list of error estimates and select the subinterval
+c           with nrmax-th largest error estimate (to be bisected next).
+c
+  110   call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax)
+c ***jump out of do-loop
+        if(errsum.le.errbnd) go to 190
+c ***jump out of do-loop
+        if(ier.ne.0) go to 170
+        if(noext) go to 160
+        erlarg = erlarg-erlast
+        if(levcur+1.le.levmax) erlarg = erlarg+erro12
+        if(extrap) go to 120
+c
+c           test whether the interval to be bisected next is the
+c           smallest interval.
+c
+        if(level(maxerr)+1.le.levmax) go to 160
+        extrap = .true.
+        nrmax = 2
+  120   if(ierro.eq.3.or.erlarg.le.ertest) go to 140
+c
+c           the smallest interval has the largest error.
+c           before bisecting decrease the sum of the errors over
+c           the larger intervals (erlarg) and perform extrapolation.
+c
+        id = nrmax
+        jupbnd = last
+        if(last.gt.(2+limit/2)) jupbnd = limit+3-last
+        do 130 k = id,jupbnd
+          maxerr = iord(nrmax)
+          errmax = elist(maxerr)
+c ***jump out of do-loop
+          if(level(maxerr)+1.le.levmax) go to 160
+          nrmax = nrmax+1
+  130   continue
+c
+c           perform extrapolation.
+c
+  140   numrl2 = numrl2+1
+        rlist2(numrl2) = area
+        if(numrl2.le.2) go to 155
+        call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres)
+        ktmin = ktmin+1
+        if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5
+        if(abseps.ge.abserr) go to 150
+        ktmin = 0
+        abserr = abseps
+        result = reseps
+        correc = erlarg
+        ertest = dmax1(epsabs,epsrel*dabs(reseps))
+c ***jump out of do-loop
+        if(abserr.lt.ertest) go to 170
+c
+c           prepare bisection of the smallest interval.
+c
+  150   if(numrl2.eq.1) noext = .true.
+        if(ier.ge.5) go to 170
+  155   maxerr = iord(1)
+        errmax = elist(maxerr)
+        nrmax = 1
+        extrap = .false.
+        levmax = levmax+1
+        erlarg = errsum
+  160 continue
+c
+c           set the final result.
+c           ---------------------
+c
+c
+  170 if(abserr.eq.oflow) go to 190
+      if((ier+ierro).eq.0) go to 180
+      if(ierro.eq.3) abserr = abserr+correc
+      if(ier.eq.0) ier = 3
+      if(result.ne.0.0d+00.and.area.ne.0.0d+00)go to 175
+      if(abserr.gt.errsum)go to 190
+      if(area.eq.0.0d+00) go to 210
+      go to 180
+  175 if(abserr/dabs(result).gt.errsum/dabs(area))go to 190
+c
+c           test on divergence.
+c
+  180 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le.
+     *  resabs*0.1d-01) go to 210
+      if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03.or.
+     *  errsum.gt.dabs(area)) ier = 6
+      go to 210
+c
+c           compute global integral sum.
+c
+  190 result = 0.0d+00
+      do 200 k = 1,last
+        result = result+rlist(k)
+  200 continue
+      abserr = errsum
+  210 if(ier.gt.2) ier = ier-1
+      result = result*sign
+  999 return
+      end