      subroutine trapzd(func,a,b,s,n) 
*
* Numerical integration by Trapezoid rule. (Called by QSIMP)
* From Numerical Recipies
*
      implicit double precision(a-h,o-z)

      external func
      logical debug
      debug = .false.

      if (n.eq.1) then
         s = 0.5*(b-a)*(func(a)+func(b))
         it = 1
       else
         tnm = it 
         del = (b-a)/tnm
         x = a + 0.5 * del
         sum = 0. 
         do 11 j = 1,it
            dum = func(x)
            sum = sum + dum 
            x = x + del
31          format(1x,2f10.5)
11       continue 
         s = 0.5 * (s + (b-a) * sum/tnm )
         it = 2 * it
        endif

        return
        end

************************************************************************

