      subroutine trapex (f,t,y,tend,eps,error,hmax,h,nstep)             
c                                                                       
c                                                                    
c     SOLUTION OF QUADRATURE PROBLEMS ON FINITE INTERVALS               
c     USING TRAPEZOIDAL SUMS WITH POLYNOMIAL EXTRAPOLATION              
c     (ROMBERG QUADRATURE)                                       
c                                                                      
c     LOCAL RELATIVE ERROR CONTROL WITH INTERNAL SCALING            
c                                                                       
c                                                                       
c     VERSION:        JUL., '82 
c     LATEST CHANGE:  JAN., '04
c                                                                       
c                                                                       
c                                                                       
c  REFERENCES:                                                          
c /1/ W.ROMBERG:                                                        
c     VEREINFACHTE NUMERISCHE INTEGRATION                               
c     DET KONG. NORSKE VIDENSKABERS SELSKABS FORHANDLINGER 28, NR.7,    
c     TRONDHEIM 1955    
c                                                                       
c /2/ R.BULIRSCH:                                                       
c     BEMERKUNGEN ZUR ROMBERG-INTEGRATION                               
c     NUM.MATH. 6, 6-16 (1964)                                          
c                                                                       
c /3/ P.DEUFLHARD:                                                      
c     ORDER AND STEPSIZE CONTROL IN EXTRAPOLATION METHODS               
c     TO APPEAR IN NUM.MATH.  (1983)                                    
c                                                                       
c /4/ P.DEUFLHARD, H.J.BAUER:                                           
c     A NOTE ON ROMBERG QUADRATURE                                      
c     UNIVERSITY OF HEIDELBERG, SFB 123: TECH. REP. 169 (1982)          
c                                                                       
c                                                                       
c                                                                       
c  EXTERNAL FUNCTION (TO BE SUPPLIED BY THE USER)                       
c                                                                       
c    F(T)               INTEGRAND                                       
c      T                ACTUAL POSITION                                 
c                                                                       
c                                                                       
c  INPUT PARAMETERS(* MARKS TRANSIENT PARAMETERS)                       
c                                                                       
c  * T                  STARTING POINT OF INTEGRATION                   
c  * Y                  INITIAL VALUE OF THE INTEGRAL                   
c                       (STANDARD INPUT VALUE : 0.)                     
c    TEND               PRESCRIBED FINAL POINT OF INTEGRATION           
c    EPS                PRESCRIBED LOCAL RELATIVE PRECISION (.GT.0)     
c    HMAX               MAXIMUM PERMITTED STEPSIZE                      
c                       (POSSIBLE STANDARD VALUE : TEND-T)              
c  * H                  INITIAL STEPSIZE PROPOSAL                       
c                       (POSSIBLE STANDARD VALUE : TEND-T)              
c                                                                       
c                                                                       
c  OUTPUT PARAMETERS                                                    
c                                                                       
c    T                  ACHIEVED FINAL POINT OF INTEGRATION             
c    Y                  NUMERICAL RESULT, IF T.EQ.TEND                  
c                        FINAL VALUE    , IF T.NE.TEND                  
c    H                  STEPSIZE PROPOSAL FOR NEXT INTEGRATION STEP     
c                       (H.EQ.0. ,IF TRAPEX FAILS TO PROCEED,           
c                        OR MORE THAN NSTMAX STEPS PER INTERVAL,        
c                        OR MORE THAN JRMAX STEPSIZE REDUCTIONS)        
c    ERROR              GLOBAL RELATIVE ERROR ESTIMATE                  
c                                                                       
c                                                                       
c  CONVERSION TO SINGLE PRECISION:                                      
c  JUST CANCEL THE CSP IN THE CSP COMMENT CARDS AND CONVERT             
c  THE ASSOCIATED PRECEDING CARDS TO CDP COMMENT CARDS                  
c                                                                       
c                                                                       
      integer nj( 7),incr( 7),nred( 6)                                  
c     INTEGER NJ(JM),INCR(JM),NRED(KM)                                  
      double precision dt( 7),d( 7, 7),a( 7),al( 7, 7)                  
cSP       REAL DT( 7),D( 7, 7),A( 7),AL( 7, 7)                          
      double precision a1,b1,c,eph,epmach,eps,err,errh,error            
cSP       REAL         A1,B1,C,EPH,EPMACH,EPS,ERR,ERRH,ERROR            
      double precision f,fc,fcm,fco,fi,fmin,g,h,half,hmax               
cSP       REAL         F,FC,FCM,FCO,FI,FMIN,G,H,HALF,HMAX               
      double precision hmaxu,hr,h1,omj,omjo,one,one1,q,red,ro           
cSP       REAL         HMAXU,HR,H1,OMJ,OMJO,ONE,ONE1,Q,RED,RO           
      double precision safe,small,t,ta,ten,tend,teps,tn,u,u1            
cSP       REAL         SAFE,SMALL,T,TA,TEN,TEND,TEPS,TN,U,U1            
      double precision v,w,y,yl,yl1,yn,yr,zero,zq,zq1                   
cSP       REAL         V,W,Y,YL,YL1,YN,YR,ZERO,ZQ,ZQ1                   
      data zero/0.d0/,fmin/1.d-2/,ro/0.25d0/,half/0.5d0/                
cSP       DATA ZERO/0.E0/,FMIN/1.E-2/,RO/0.25E0/,HALF/0.5E0/            
      data safe/0.5d0/,one/1.d0/,one1/1.01d0/,ten/1.d1/                 
cSP       DATA SAFE/0.5E0/,ONE/1.E0/,ONE1/1.01E0/,TEN/1.E1/             
      data dt/ 7*0.d0/                                                  
cSP       DATA DT/ 7*0.E0/                                              
c     DATA DT/JM*0./                                                    
c                                                                       
c  STEPSIZE SEQUENCE  ( DUE TO /1/ )                                    
      data nj/1,2,3,4,6,8,12/                                           
c  ASSOCIATED NORMALIZED WORK PER UNIT STEP                             
      data a/1.d0,2.d0,4.d0,6.d0,8.d0,12.d0,16.d0/                      
cSP       DATA A/1.E0,2.E0,4.E0,6.E0,8.E0,12.E0,16.E0/                  
c                                                                       
      call zibconst(epmach, small)
c                                                                       
c  MAXIMUM COLUMN NUMBER (1.LE.KM.LE.6)                                 
      km=6                                                              
c                                                                       
c  ASSOCIATED MAXIMUM ROW NUMBER (2.LE.JM.LE.7)                         
      jm=km+1                                                           
c                                                                       
c  INTERNAL PARAMETERS                                                  
c  STANDARD VALUES FIXED BELOW                                          
c                                                                       
c                                                                       
c  MAXIMUM PERMITTED NUMBER OF INTEGRATION STEPS PER INTERVAL           
      nstmax=10000                                                      
c                                                                       
c  MAXIMUM PERMITTED NUMBER OF STEPSIZE REDUCTIONS                      
      jrmax=50                                                          
c                                                                       
c  INITIAL PREPARATIONS                                                 
      eph=ro*eps                                                        
      error=zero                                                        
      omjo=zero                                                         
      a1=a(1)                                                           
      u1=zero                                                           
      do 12 j=2,jm                                                      
      j1=j-1                                                            
      incr(j1)=0                                                        
      nred(j1)=0                                                        
      v=dfloat(nj(j))                                                   
cSP       V=FLOAT(NJ(J))                                                
      do 10 k=1,j1                                                      
      w=v/dfloat(nj(k))                                                 
cSP       W=V/FLOAT(NJ(K))                                              
10    d(j,k)=w*w                                                        
      if(j.eq.2) goto 12                                                
      v=a(j)                                                            
      w=v-a1+one                                                        
      do 11 k1=2,j1                                                     
      k=k1-1                                                            
      u=(a(k1)-v)/(w*dfloat(k+k1))                                      
cSP       U=(A(K1)-V)/(W*FLOAT(K+K1))                                   
      u=eph**u
      al(j1,k)=u                                                        
      u=u*a(k1)                                                         
      if(k1.gt.2.and.u*one1.gt.u1) goto 13                              
11    u1=u                                                              
      if(a(j)*one1.gt.u1) goto 13                                       
12    continue                                                          
      koh=km                                                            
      joh=jm                                                            
      goto 14                                                           
13    joh=max0(2,j1-1)                                                  
      koh=joh-1                                                         
14    km=koh                                                            
      jm=km+1                                                           
      incr(jm)=-1                                                       
      epmach=epmach*ten                                                 
      hmax=dabs(hmax)                                                   
cSP       HMAX=ABS(HMAX)                                                
      nstep=0                                                           
      teps=(dabs(t)+dabs(tend))*epmach                                  
cSP       TEPS=(ABS(T)+ABS(TEND))*EPMACH                                
      h1=tend-t                                                         
      if(dabs(h1).le.teps) return                                       
cSP       IF(ABS(H1).LE.TEPS) RETURN                                    
      q=h1/h                                                            
      if(q.le.epmach) return                                            
      hmaxu=hmax                                                        
      hr=hmax                                                           
      zq=f(t)                                                           
      yr=zero                                                           
c                                                                       
15    if(q.ge.one1) goto 16                                             
      hr=h                                                              
      h=h1                                                              
16    jred=0                                                            
      do 17 k=1,km                                                      
17    incr(k)=incr(k)+1                                                 
      hmax=dabs(h1)                                                     
cSP       HMAX=ABS(H1)                                                  
      if(hmaxu.lt.hmax) hmax=hmaxu                                      
c                                                                       
c  TRAPEZOIDAL SUM                                                      
c                                                                       
20    tn=t+h                                                            
      if(h1.eq.h) tn=tend                                               
      fcm=dabs(h)/hmax                                                  
cSP       FCM=ABS(H)/HMAX                                               
      if(fcm.lt.fmin) fcm=fmin                                          
c                                                                       
      yn=zero                                                           
      yl=zero                                                           
      do 35 j=1,jm                                                      
      m=nj(j)                                                           
      g=h/dfloat(m)                                                     
cSP       G=H/FLOAT(M)                                                  
      if(j.gt.1) goto 21                                                
      zq1=f(t+h)                                                        
      yl=(zq1+zq)*half                                                  
      yr=yl*g                                                           
      goto 30                                                           
21    if(mod(m,3).eq.0) goto 23                                         
      yl1=yl                                                            
      do 22 i=1,m,2                                                     
      yl=yl+f(t+dfloat(i)*g)                                            
cSP       YL=YL+F(T+FLOAT(I)*G)                                         
22    continue                                                          
      yr=yl*g                                                           
      goto 30                                                           
23    do 24 i=1,m,6                                                     
      fi=dfloat(i)                                                      
cSP       FI=FLOAT(I)                                                   
24    yn=yn+f(t+fi*g)+f(tn-fi*g)                                        
      yr=(yn+yl1)*g                                                     
c                                                                       
c  EXTRAPOLATION                                                        
c                                                                       
30    v=dt(1)                                                           
      c=yr                                                              
      dt(1)=c                                                           
      if(j.eq.1) goto 35                                                
      ta=c                                                              
      do 31 k=2,j                                                       
      jk=j-k+1                                                          
      b1=d(j,jk)                                                        
      w=c-v                                                             
      u=w/(b1-one)                                                      
      c=b1*u                                                            
      v=dt(k)                                                           
      dt(k)=u                                                           
31    ta=u+ta                                                           
      errh=dabs(u)                                                      
cSP       ERRH=ABS(U)                                                   
      yr=ta                                                             
      ta=dabs(y+yr)                                                     
cSP       TA=ABS(Y+YR)                                                  
      if(ta.lt.small) ta=one                                            
      err=errh/ta                                                       
      konv=0                                                            
      if(err.lt.eps) konv=1                                             
      err=err/eph                                                       
c                                                                       
c ORDER CONTROL                                                         
      k=j-1                                                             
      l=j+k                                                             
      fc=err**(one/dfloat(l))                                           
cSP       FC=ERR**(ONE/FLOAT(L))                                        
      if(fc.lt.fcm) fc=fcm                                              
c  OPTIMAL ORDER DETERMINATION                                          
      omj=fc*a(j)                                                       
      if(j.gt.2.and.omj*one1.gt.omjo.or.k.gt.joh) goto 32               
      ko=k                                                              
      jo=j                                                              
      omjo=omj                                                          
      fco=fc                                                            
32    if(j.lt.koh.and.nstep.gt.0) goto 35                               
      if(konv.eq.0) goto 33                                             
      if(ko.lt.k.or.incr(j).lt.0) goto 40                               
c  POSSIBLE INCREASE OF ORDER                                           
      if(nred(ko).gt.0) nred(ko)=nred(ko)-1                             
      if(j.ge.jm) goto 40                                               
      fc=fco/al(j,k)                                                    
      if(fc.lt.fcm) fc=fcm                                              
      j1=j+1                                                            
      if(a(j1)*fc*one1.gt.omjo) goto 40                                 
      fco=fc                                                            
      ko=jo                                                             
      jo=jo+1                                                           
      goto 40                                                           
c                                                                       
c                                                                       
c  CONVERGENCE MONITOR                                                  
33    red=one/fco                                                       
      jk=km                                                             
      if(joh.lt.km) jk=joh                                              
      if(k.ge.jk) goto 36                                               
      if(ko.lt.koh) red=al(koh,ko)/fco                                  
34    if(al(jk,ko).lt.fco) goto 36                                      
35    continue                                                          
c                                                                       
c STEPSIZE REDUCTION                                                    
36    red=red*safe                                                      
      h=h*red                                                           
      if(nstep.eq.0) goto 38                                            
      nred(koh)=nred(koh)+1                                             
      do 37 l=koh,km                                                    
37    incr(l)=-2-nred(koh)                                              
38    jred=jred+1                                                       
      if(jred.gt.jrmax) goto 51                                         
      goto 20                                                           
c                                                                       
c  PREPARATIONS FOR NEXT BASIC STEP                                     
40    t=tn
      h1=tend-t                                                         
      y=y+yr                                                            
      error=error+errh                                                  
      zq=zq1                                                            
      nstep=nstep+1
c      write(6,*) nstep,t,y,error,ko
      if(nstep.gt.nstmax) go to 51                                      
c                                                                       
c STEPSIZE PREDICTION                                                   
      if(fco.ne.fcm) hr=h                                               
      h=h/fco                                                           
      koh=ko                                                            
      joh=koh+1                                                         
      if(dabs(h).le.teps) go to 51                                      
cSP       IF(ABS(H).LE.TEPS) GO TO 51                                   
c                                                                       
41    continue                                                          
      if(dabs(h1).le.teps) goto 50                                      
cSP       IF(ABS(H1).LE.TEPS) GOTO 50                                   
      q=h1/h                                                            
      if(q.gt.epmach) goto 15                                           
c                                                                       
c  SOLUTION EXIT                                                        
50    h=hr                                                              
      ta=dabs(y)                                                        
cSP       TA=ABS(Y)                                                     
      if(ta.lt.small) ta=small                                          
      error=error/ta                                                    
      hmax=hmaxu                                                        
      return                                                            
c                                                                       
c  FAIL EXIT                                                            
51    h=zero                                                            
      hmax=hmaxu                                                        
      return                                                            
c                                                                       
c  END TRAPEX                                                           
c                                                                       
      end                                                               
