program romb1
c
c -- numerical integration by the romberg method
c
      integer out
      real sum, upper, lower, tol
      external f
      common /inout/ out
      data lower/1.0/, upper/9.0/, tol/1.0e-5/
c
      out = 6
      write(out,*) ' romberg integration'
      call romber(lower, upper, tol, sum, f)
      write(out, 104) sum
      stop
104   format(/' area =', f10.5/)
      end
      function f(x)
c
      f = 1.0 / x
      return
      end
      subroutine romber(lower, upper, tol, ans, f)
c
c -- numerical integration by the romberg method
c
      integer out, pieces, i, nx(16)
      integer nt, ii, n, nn, l, ntra, k, m, j, l2
      real x, delta, lower, upper, sum, tol
      real c, fotom, t(136), ans
      common /inout/ out
c
      pieces  = 1
      nx(1) = 1
      delta = (upper - lower) / pieces
      c = (f(lower) + f(upper)) / 2.0
      sum = c
      t(1) = delta * c
      n = 1
      nn = 2
5     n = n + 1
       fotom = 4.0
       nx(n) = nn
       pieces = pieces * 2
       l = pieces - 1
       l2 = (l+1) / 2
       delta = (upper - lower) / pieces
c -- compute trapezoidal sum for 2**(n-1)+1 points
       do 10 ii = 1, l2
        i = ii * 2 - 1
        x = lower + delta * i
        sum = sum + f(x)
10     continue
       t(nn) = sum * delta
       ntra = nx(n-1)
       k = n-1
c -- compute n-th row ot t array
       do 20 m = 1, k
        j = nn + m
        nt = nx(n-1) + m - 1
        t(j) = (fotom * t(j-1) - t(nt)) / (fotom - 1)
        fotom = fotom * 4
20     continue
       write(out, 101) pieces, t(nn), j, t(j)
       if (n .lt. 5)  goto 40
        if (t(nn+1) .eq. 0.0) goto 40
         if (abs(t(ntra+1)-t(nn+1)) .le. abs(t(nn+1)*tol)) goto 99
         if (abs(t(nn-1) - t(j)) .le. abs(t(j)*tol)) goto 99
         if (n .gt. 15) goto 100
40     nn = j + 1
      goto 5
99    ans = t(j)
      return
100   write(out,*) ' romberg error'
      return
101   format(1x, i7, f9.5, i5, f9.5)
      end