       subroutine qromb ( func, x1, x2, ss )
c
c    dac:zeus3d.qromb <------------------- integrates func from x1 to x2
c    from whp:numerical recipes                           december, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  Returns the definite integral of the function "func"
c  between specified limits "x1" and "x2" using Romberg's method of
c  order 2k, where k=2 is Simpson's rule (see Numerical Recipes, 1st
c  edition for FORTRAN, page 114).
c
c  INPUT VARIABLES:
c    func     name of external function describing the univariate
c             function to be integrated.
c    x1, x2   integration limits
c
c  OUTPUT VARIABLES:
c    ss       value of definite integral
c
c  EXTERNALS:
c    FUNC
c    TRAPZD
c    POLINT
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c------------------------ IMPLICIT STATEMENTS --------------------------
c-----------------------------------------------------------------------
c
       implicit      none
c
c-----------------------------------------------------------------------
c----------------------------- PARAMETERS ------------------------------
c-----------------------------------------------------------------------
c
c      Primary parameters:
c
       integer       in      , jn      , kn      , ijkn    , nxpx
     1             , nypx    , nxrd    , nyrd    , niov    , nios
     2             , ntsl    , mn      , nmd     , isig
       real*8        pi  , tiny  , huge  , one  , mone  , haf  , zro
c
       parameter     ( in = 5005, jn = 1, kn = 10, ijkn = 5005)
       parameter     ( nxpx = 1, nypx = 1, nxrd =   1, nyrd =   1 )
       parameter     ( niov =  20, nios =  20, ntsl =1000 )
       parameter     ( mn   =   1, nmd  =   1 )
       parameter     ( pi   = 3.1415926535898 )
       parameter     ( isig = 9, tiny = 1.0d-99, huge = 1.0d+99 )
       parameter     ( one  = 1.0d+00, mone = -1.0d+00
     1                               , haf  =  0.5d+00, zro = 0.0d+00)
c
c      Secondary parameters:
c
       integer        idim    , jdim
     1              , nt01    , nt02    , nt03    , nt04    , nt05
     2              , nt06    , nt07    , nt08    , nt09    , nt10
     3              , nt11    , nt12    , nt13    , nt14    , nt15
     4              , nt16    , nt17    , nt18    , nt19    , nt20
       parameter     ( idim =   kn, jdim =   in )
       parameter     ( nt01 = 600, nt02 = 600, nt03 =   7, nt04 =   5 )
       parameter     ( nt05 =  18, nt06 =   1, nt07 =  10, nt08 =  12 )
       parameter     ( nt09 =   5, nt10 =  54, nt11 =   8 )
       parameter     ( nt12 =        16, nt17 =         4 )
       parameter     ( nt13 = nt12 +  4, nt18 = nt17 +  1 )
       parameter     ( nt14 = nt13     , nt19 = nt18      )
       parameter     ( nt15 = nt14 + 30, nt20 = nt19 +  3 )
       parameter     ( nt16 = nt14 + 42 )
c
       integer       jmx     , jmxp    , k       , km
       real*8        eps
       parameter     ( jmx=20, jmxp=jmx+1, k=5, km=k-1 )
       parameter     ( eps=0.1**(isig+1) )
c
       integer       j
       real*8        x1      , x2      , ss      , dss
c
       real*8        s       (jmxp), h       (jmxp)
c
c      External statements
c
       real*8        func
       external      func    , trapzd  , polint
c
c-----------------------------------------------------------------------
c
       h(1) = 1.0
       do 10 j=1,jmx
         call trapzd ( func, x1, x2, j, s(j) )
         if (j .ge. k) then
           call polint ( h(j-km), s(j-km), k, zro, ss, dss )
           if (abs(dss) .lt. eps*abs(ss)) return
         endif
         s(j+1) = s(j)
         h(j+1) = 0.25 * h(j)
10     continue
       write (6, 2000) jmx, dss, ss
       return
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2000   format('QROMB   : *** WARNING *** Romberg Integration failed '
     1       ,'to converge in ',i2,' steps.',/
     1       ,'QROMB   : error =',1pg12.5,', definite integral =',g12.5)
c
       end
