From d21f9b12f47ca982384b99aabab09e56cb18f4b7 Mon Sep 17 00:00:00 2001 From: Daniela Farina Date: Wed, 27 Nov 2013 14:27:02 +0000 Subject: [PATCH] print in fort.33 also for tau>taucr; changed test for Npl>1 in after_onestep --- Makefile | 2 +- src/gray.f | 10 +- src/grayl.f | 777 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 783 insertions(+), 6 deletions(-) diff --git a/Makefile b/Makefile index 173e322..160f0ae 100644 --- a/Makefile +++ b/Makefile @@ -11,7 +11,7 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -FFLAGS=-O3 +FFLAGS=-O3 DIRECTIVES = -DREVISION="'$(shell svnversion src)'" diff --git a/src/gray.f b/src/gray.f index 0805c49..6f368da 100644 --- a/src/gray.f +++ b/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 diff --git a/src/grayl.f b/src/grayl.f index d1aa7b9..0f0a9f6 100644 --- a/src/grayl.f +++ b/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