subroutine romint ( val, err, eps, a, b, n, maxe, f )

c*********************************************************************72
c
cc romint
c
      implicit none

      real a
      real b
      real bb
      real eps
      real err
      external f
      real f
      real h
      integer j
      integer k
      integer k0
      integer k1
      integer k2
      integer kk
      integer kkk
      integer l
      integer maxe
      integer n
      integer n0
      integer n1
      real r
      real rm(16)
      real s
      real s0
      real s1
      real t
      real val
c
c  initial trapezoid rule.
c
      t = ( b - a ) * ( f(a) + f(b) ) * 0.5e+00
c
c  initial rectangle value.
c
      rm(1) = ( b - a ) * f ( ( a + b ) * 0.5e+00 )

      n = 2
      r = 4.0e+00

      do 11 k = 1, maxe

        bb = ( r * 0.5e+00 - 1.0e+00 ) / ( r - 1.0e+00 )
c
c  improved trapezoid value.
c
        t = rm(1) + bb * ( t - rm(1) )
c
c  double number of subdivisions of (a,b).
c
        n = 2 * n
        s = 0.0e+00
        h = ( b - a ) / float ( n )
c
c  calculate rectangle value.
c
        if ( n - 32 ) 1, 1, 2
1       n0 = n
        go to 4
2       n0 = 32
        if ( n - 512 ) 4, 4, 5
4       n1 = n
        go to 6
5       n1 = 512
6       do 9 k2 = 1, n, 512
          s1 = 0.0e+00
          kk = k2 + n1 - 1
          do 8 k1 = k2, kk, 32
            s0 = 0.0e+00
            kkk = k1 + n0 - 1
            do k0 = k1, kkk, 2
              s0 = s0 + f ( a + float ( k0 ) * h )
            end do
            s1 = s1 + s0
8         continue
          s = s + s1
9       continue
        rm(k+1) = 2.0e+00 * h * s
c
c  end calculation of rectangle value.
c
        r = 4.0e+00
c
c  form romberg table from rectangle values.
c
        do j = 1, k
          l = k + 1 - j
          rm(l) = rm(l+1) + ( rm(l+1) - rm(l) ) / ( r - 1.0e+00 )
          r = 4.0e+00 * r
        end do

        err = abs ( t - rm(1) ) * 0.5e+00
c
c  convergence test.
c
        if ( err - eps ) 12, 12, 11

11    continue

      k = 0

12    val = ( t + rm(1) ) * 0.5e+00
      n = n + 1
      maxe = k

      return
      end