program trap1
c
      integer out, pieces
      real sum, upper, lower
      data lower/1.0/, upper/9.0/
c
      out = 6
      write(out, *) ' trapezoidal integration'
10    write(out, 102)
      read(*, 103) pieces
      if (pieces .lt. 0) goto 99
      call trapez(lower, upper, pieces, sum)
      write(out, 104) sum
      goto 10
99    stop
102   format(' how many pieces? ' )
103   format(i3)
104   format(' area =', f10.5)
      end
      subroutine trapez(lower, upper, pieces, sum)
c
c -- numerical integration by the trapezoid method
c
      integer pieces, i
      real x, delta, endsum, midsum, lower, upper, sum
c
c -- f(x) = 1 / x, be careful of x = 0
c
      f(x) = 1.0 / x
c
      delta = (upper - lower) / pieces
      endsum = f(lower) + f(upper)
      midsum = 0.0
      do 10 i = 1, pieces
       x = lower + i * delta
       midsum = midsum + f(x)
10    continue
      sum = (endsum + 2.0*midsum) * delta * 0.5
      return
      end