       subroutine trapzd ( func, x1, x2, n, s )
c
c    dac:zeus3d.trapzd <-------------- n'th refinement of trapezoid rule
c    from whp:numerical recipes                           december, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  This routine returns in "s" the n'th stage of refinement of
c  an extended trapezoidal rule.  "func" is input as the name of the
c  function to be integrated between limits "x1" and "x2", also input.
c  When called with n=1, the crudest estimate of the definite integral
c  is returned.  Subsequent calls with n=2,3,... (in that order) will
c  improve the accuracy of "s" by adding 2**(n-2) additional interior
c  points.  The value of "s" should not be modified between successive
c  calls (see Numerical Recipes, 1st edition for FORTRAN, page 111).
c
c  INPUT VARIABLES:
c    func     name of external function describing the univariate
c             function to be integrated.
c    x1, x2   integration limits
c    n        stage of refinement
c
c  OUTPUT VARIABLES:
c    s        value of definite integral after n'th stage of refinement.
c
c  EXTERNALS: [NONE]
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c------------------------ IMPLICIT STATEMENTS --------------------------
c-----------------------------------------------------------------------
c
       implicit      none
c
       integer       n       , it      , j
       real*8        x1      , x2      , s      , tnm     , del
     1             , x       , sum
c
c      Save statements
c
       save          it
c
c      External statements
c
       real*8        func
       external      func
c
c-----------------------------------------------------------------------
c
       if (n .eq. 1) then
         s  = 0.5 * ( x2 - x1 ) * ( func(x2) + func(x1) )
         it = 1
       else
         tnm = 1.0 /  real(it)
         del = ( x2 - x1 ) * tnm
         x   = x1 + 0.5 * del
         sum = 0.0
         do 10 j=1,it
           sum = sum + func(x)
           x   = x + del
10       continue
         s  = 0.5 * ( s + ( x2 - x1 ) * sum * tnm )
         it = 2 * it
       endif
c
       return
       end
