! aquad3e.f90 from book page 301, with modifications 

module global
  integer, parameter :: maxtop = 100
  double precision, dimension(maxtop,7) :: stack  ! a place to store/retrieve 
  integer :: top                                  ! top of stack 
end module global

subroutine store(s1, s2, s3, s4, s5, s6, s7)
  use global
  implicit none
  double precision, intent(in) :: s1, s2, s3, s4, s5, s6, s7
  stack(top,1) = s1
  stack(top,2) = s2
  stack(top,3) = s3
  stack(top,4) = s4
  stack(top,5) = s5
  stack(top,6) = s6
  stack(top,7) = s7
end subroutine store

subroutine retrieve(s1, s2, s3, s4, s5, s6, s7)
  use global
  implicit none
  double precision, intent(out) :: s1, s2, s3, s4, s5, s6, s7
  s1 = stack(top,1)
  s2 = stack(top,2)
  s3 = stack(top,3)
  s4 = stack(top,4)
  s5 = stack(top,5)
  s6 = stack(top,6)
  s7 = stack(top,7)
end subroutine retrieve 

function Sn(F0, F1, F2, h) result(integral) ! 2h/6
  implicit none
  double precision, intent(in) :: F0, F1, F2, h
  double precision integral
  integral = h*(F0 + 4.0*F1 + F2)/3.0
  ! error term 1/12 h^3 f^(2)(c) 
end function Sn 

function RS(F0, F1, F2, F3, F4, h) result(integral) ! 180h/45 
  implicit none
  double precision, intent(in) :: F0, F1, F2, F3, F4, h
  double precision integral
  integral = h*(14.0*F0 +64.0*F1 + 24.0*F2 + 64.0*F3 + 14.0*F4)/45.0 
  ! error term  8/945  h^7 f^(8)(c) 
end function RS 

function aquad3e(f, xmin, xmax, eps, i, j) result(integral)
  use global
  implicit none
  double precision, intent(in) :: xmin, xmax, eps
  double precision integral
  integer, intent(in) :: i, j
  double precision :: a, b, c, d, e
  double precision :: Fa, Fb, Fc, Fd, Fe
  double precision :: h1, h2, hs
  double precision :: Sab, Sac, Scb, S2ab
  double precision :: tol ! eps 
  double precision :: val, value
  integer, parameter :: ns = 9
  integer :: k, nf

  interface
    function f(x,i,j) result(y)
      double precision, intent(in) :: x
      integer, intent(in) :: i, j
      double precision :: y
    end function f
    subroutine store(s1, s2, s3, s4, s5, s6, s7)
      double precision, intent(in) :: s1, s2, s3, s4, s5, s6, s7
    end subroutine store
    subroutine retrieve(s1, s2, s3, s4, s5, s6, s7)
      double precision, intent(out) :: s1, s2, s3, s4, s5, s6, s7
    end subroutine retrieve
    function Sn(F0, F1, F2, h) result(integral)
      double precision, intent(in) :: F0, F1, F2, h
      double precision integral
    end function Sn 
    function RS(F0, F1, F2, F3, F4, h) result(integral) 
      double precision, intent(in) :: F0, F1, F2, F3, F4, h
      double precision integral
    end function RS 
  end interface

  ! print *, 'aquad3f i=',i,', j=',j
  nf = 0
  top = 1
  value = 0.0
  tol = eps
  hs = (xmax-xmin)/ns
  a = xmin
  b = a + hs
  do k=1,ns
    h1 = (b-a)/2.0
    c = a + h1
    Fa = f(a,i,j)
    Fc = f(c,i,j)
    Fb = f(b,i,j)
    nf = nf + 3
    Sab = Sn(Fa, Fc, Fb, h1)
    call store(a, Fa, Fc, Fb, h1, tol, Sab)
    top = top+1
    a = b
    b = a + hs
  end do
  do while(top .gt. 1)
    top = top-1
    call retrieve(a, Fa, Fc, Fb, h1, tol, Sab)
    c = a + h1
    b = a + 2.0*h1
    h2 = h1/2
    d = a + h2
    e = a + 3.0*h2
    Fd = f(d,i,j)
    Fe = f(e,i,j)
    nf = nf + 2
    Sac = Sn(Fa, Fd, Fc, h2)
    Scb = Sn(Fc, Fe, Fb, h2)
    S2ab = Sac + Scb
    if(abs(S2ab-Sab) .lt. tol) then
      val = RS(Fa, Fd, Fc, Fe, Fb, h2)
      value = value + val
    else
      h1 = h2
      tol = tol/2.0
      call store(a, Fa, Fd, Fc, h1, tol, Sac)
      top = top+1
      call store(c, Fc, Fe, Fb, h1, tol, Scb)
      top = top+1
    end if
    if(top .ge. maxtop) exit
  end do 
  integral = value
  ! print *, 'aquad3e finished, number of function evaluation=', nf
end function aquad3e