cMM this is the original adapt code with one modification.
cMM instead of calling the external function "functn", a fixed
cMM external routine adphlp is always called, and passed a pointer
cMM to the external s function.
cMM					Michael Meyer, October 1989.

      subroutine adapt(ndim,a,b,minpts,maxpts,eps,relerr,
     *	   lenwrk,wrkstr,finest,ifail)
c***begin prologue adapt
c  adaptive multidimensional integration subroutine
c	    author: A. C. Genz, Washington State University
c		     19 March 1984
c**************	 parameters for adapt  ********************************
c***** input parameters
c  ndim	   number of variables, must exceed 1, but not exceed 20
c  a	   real array of lower limits, with dimension ndim
c  b	   real array of upper limits, with dimension ndim
c  minpts  minimum number of function evaluations to be allowed.
c	   on the first call to adapt minpts should be set to a
c	   non negative value (caution... minpts is altered by adapt).
c	   It is possible to continue a calculation to greater accuracy
c	   by calling adapt again by decreasing eps (described below)
c	   and resetting minpts to any negative value.
c	   minpts must not exceed maxpts.
c  maxpts  maximum number of function evaluations to be allowed,
c	   which must be at least rulcls, where
c	   rulcls =  2**ndim+2*ndim**2+6*ndim+1
c
c	       for ndim =  2   3   4   5   6   7   8   9   10	12    15      20
c      maxpts >= rulcls = 25  45  73 113 173 269 433 729 1285 4457 33309 1049497
c
c	  a SUGGESTED value for maxpts is 100 times the above values.
c
c  functn  externally declared user defined function to be integrated.
c	   it must have parameters (ndim,z), where z is a real array
c	   of dimension ndim.
cTSL this function has been replaced by the fixed function adhlp
c  eps	   required relative accuracy
c  lenwrk  length of array wrkstr of working storage, the routine
c	   needs (2*ndim+3)*(1+maxpts/rulcls)/2 for lenwrk if
c	   maxpts function calls are used.
c	   for guidance, if you set maxpts to 100*rulcls (see table
c	   above) then acceptable values for lenwrk are
c
c	     for ndim = 2    3	  4    5    6	 7    8	    9
c	     lenwrk =  357  561	 1785 3417 6681 13209 26265 52377
c
c***** OUTPUT parameters
c
c  minpts  actual number of function evaluations used by adapt
c  wrkstr  real array of working storage of dimension (lenwrk).
c  relerr  estimated relative accuracy of finest
c  finest  estimated value of integral ["FINal ESTimate"]
c  ifail : return code
c
c	ifail=0 for normal exit, when estimated relative accuracy
c		relerr is less than eps with maxpts or less function
c		calls made.
c	ifail=1 if maxpts was too small for adapt to obtain the
c		required relative accuracy eps.
c		In this case adapt returns a value of finest
c		with estimated relative accuracy relerr.
c	ifail=2 if lenwrk too small for maxpts function calls.
c		In this case adapt returns a value of finest with
c		estimated accuracy relerr using the working storage
c		available, but relerr will be greater than eps.
c	ifail=3 if ndim < 2, ndim > 20,
c	ifail=4 if minpts > maxpts,
c	ifail=5 if maxpts < rulcls or other memory problems
c		      (which will only be found later)
c***********************************************************************
c***end prologue adapt

      implicit none

C-- Arguments:
C      double precision functn
C      external functn

      integer ndim, minpts,maxpts, lenwrk, ifail
      double precision a(ndim), b(ndim), eps, relerr, wrkstr(lenwrk),
     &	   finest

C-- Local Variables:
      double precision center(20), width(20)
      double precision errmin, rgnerr, rgnval, half, zero,one,two

      integer divaxo, divaxn, divflg, funcls, index1, index2,
     *     j, k,  maxcls, rgnstr, rulcls, sbrgns, sbtmpp, subrgn, subtmp

      data zero/0d0/, one/1d0/, two/2d0/

c Check arguments; fail w/ code '3' or '4'
      relerr=one
      funcls=0
      ifail=3
      if(ndim.lt.2 .or. ndim.gt.20) goto 990
      ifail=4
      if(minpts.gt.maxpts) goto 990
      ifail=5
c
c*****	initialisation of subroutine
c
      half=one/two
      rgnstr = 2*ndim+3
      errmin = zero
      maxcls = 2**ndim + 2*ndim**2 + 6*ndim+1
      maxcls = min0(maxcls,maxpts)
      divaxo=0
c
c*****	end subroutine initialisation

      if(minpts.lt.0) then
c	on re-entry {with intact  lenwrk, wrkstr(1:lenwrk)} :
         sbrgns=wrkstr(lenwrk-1)
         goto 280
      endif

      do 30 j=1,ndim
	 width(j)=(b(j)-a(j))*half
 30	 center(j)=a(j)+width(j)
      finest=zero
      wrkstr(lenwrk)=zero
      divflg=1
      subrgn=rgnstr
      sbrgns=rgnstr

C-- REPEAT --- (outermost loop) -------
   40 call bsrl(ndim,center,width,maxcls,rulcls,
     *	   errmin,rgnerr,rgnval,divaxo,divaxn)
      finest=finest+rgnval
      wrkstr(lenwrk)=wrkstr(lenwrk)+rgnerr
      funcls = funcls + rulcls
c
c*****	place results of basic rule into partially ordered list
c*****	according to subregion error

      if(divflg .eq. 0) then
c
c   ***	when divflg=0 start at top of list and move down list tree to
c       find correct position for results from first half of recently
c       divided subregion
 200	 subtmp=2*subrgn
	 if(subtmp.le.sbrgns) then
	    if(subtmp.ne.sbrgns) then
	       sbtmpp=subtmp+rgnstr
	       if(wrkstr(subtmp).lt.wrkstr(sbtmpp)) subtmp=sbtmpp
	    endif
	    if(rgnerr.lt.wrkstr(subtmp)) then
	       do 220 k=1,rgnstr
		  index1=subrgn-k+1
		  index2=subtmp-k+1
		  wrkstr(index1)=wrkstr(index2)
 220	       continue
	       subrgn=subtmp

	       goto 200
	    endif
	 endif

      else
c
c   *** when divflg=1 start at bottom right branch and move up list
c       tree to find correct position for results from second half of
c       recently divided subregion
 230	 subtmp=(subrgn/(rgnstr*2))*rgnstr
	 if(subtmp.ge.rgnstr) then
	    if(rgnerr.gt.wrkstr(subtmp)) then
	       do 240 k=1,rgnstr
		  index1=subrgn-k+1
		  index2=subtmp-k+1
		  wrkstr(index1)=wrkstr(index2)
 240	       continue
	       subrgn=subtmp
	       goto 230
	    endif
	 endif
      endif

c**** store results (length= 2*ndim+3) of basic rule
c     in correct position in list:
  250 wrkstr(subrgn)=rgnerr
      wrkstr(subrgn-1)=rgnval
      wrkstr(subrgn-2)=divaxn
      do 260 j=1,ndim
	subtmp=subrgn-2*(j+1)
	wrkstr(subtmp+1)=center(j)
	wrkstr(subtmp)=width(j)
 260  continue

      if(divflg .eq. 0) then
c***  when divflg=0 prepare for second application of basic rule
	 center(divaxo)=center(divaxo)+two*width(divaxo)
	 sbrgns=sbrgns+rgnstr
	 subrgn=sbrgns
	 divflg=1
c***  loop back to apply basic rule to other half of subregion
	 go to 40
      endif
c
c*****	end ordering and storage of basic rule results
c*****	make checks for possible termination of routine
c
  270 relerr=one
      if(wrkstr(lenwrk).le.zero) wrkstr(lenwrk)=zero
      if(dabs(finest).ne.zero) relerr=wrkstr(lenwrk)/dabs(finest)
      if(relerr.gt.one) relerr=one

      if(sbrgns+rgnstr .gt. lenwrk-2) ifail=2
c     MM: Note that rgnstr = 2*ndim+3, and lenwrk
c     --  lenwrk = (2*ndim + 3)* (1 + maxpts/rulcls)/2
c         where rulcls = 2^ndim + 2*ndim^2 + 6*ndim + 1
c         and   maxpts increases (for each .C() call from R)
      if(funcls+funcls*rgnstr/sbrgns.gt.maxpts) ifail=1
      if(relerr.lt.eps.and.funcls.ge.minpts) ifail=0
      if(ifail.lt.3) goto 990
c
c*****	prepare to use basic rule on each half of subregion with largest
c	error
  280 divflg=0
      subrgn=rgnstr
      subtmp = 2*sbrgns/rgnstr
      maxcls = maxpts/subtmp
      errmin = dabs(finest)*eps/dfloat(subtmp)
      wrkstr(lenwrk)=wrkstr(lenwrk)-wrkstr(subrgn)
      finest=finest-wrkstr(subrgn-1)
      divaxo=wrkstr(subrgn-2)
      do 290 j=1,ndim
	subtmp=subrgn-2*(j+1)
	center(j)=wrkstr(subtmp+1)
  290	width(j)=wrkstr(subtmp)
      width(divaxo)=width(divaxo)*half
      center(divaxo)=center(divaxo)-width(divaxo)
c
c*****	loop back to apply basic rule
c
      goto 40
c
c*****	termination point
c
  990 minpts=funcls
      wrkstr(lenwrk-1)=sbrgns
      return
      end