File: /nfs/home/0/users/jenkins/mfix.git/model/dqmom/odeint.f

1     
2            SUBROUTINE odeint(ystart1,nvar,x1,x2,eps,h1,hmin,nok,nbad)
3            IMPLICIT NONE
4     
5            INTEGER nbad,nok,nvar,KMAXX,MAXSTP,NMAX
6     
7            DOUBLE PRECISION  eps,h1,hmin,x1,x2,ystart1(nvar),TINY
8            PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=200,TINY=1.e-30)
9            INTEGER i,kmax1,kount,nstp
10            DOUBLE PRECISION dxsav,h,hdid,hnext,x,xsav,dydx(NMAX),&
11            xp(KMAXX),y(NMAX),yp(NMAX,KMAXX),yscal(NMAX)
12            x=x1
13            h=sign(h1,x2-x1)
14            nok=0
15            nbad=0
16            kount=0
17     !#######################################
18     !      Don't store the intermediate results.
19     !#######################################
20            kmax1=0.0
21            dxsav=0.0
22     !##############################################
23     ! store the intermediate results. uncooment them.
24     !      kmax1=KMAXX
25     !       dxsav=0.0
26     !##################################################
27            do i=1,nvar
28             y(i)=ystart1(i)
29            end do
30     
31     
32            if (kmax1>0) xsav=x-2.*dxsav
33            do  nstp=1,MAXSTP
34               call source_population_eq(x,y,dydx)
35     
36              do i=1,nvar
37               yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
38              end do
39              if(kmax1>0)then
40               if(abs(x-xsav)>abs(dxsav)) then
41                 if(kount<kmax1-1)then
42                   kount=kount+1
43                   xp(kount)=x
44                   do i=1,nvar
45                     yp(i,kount)=y(i)
46                   end do
47                   xsav=x
48                  endif
49                endif
50               endif
51              if((x+h-x2)*(x+h-x1)>0.) h=x2-x
52              call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext)
53              if(hdid==h)then
54                nok=nok+1
55              else
56               nbad=nbad+1
57              endif
58              if((x-x2)*(x2-x1)>=0.)then
59               do i=1,nvar
60                 ystart1(i)=y(i)
61               end do
62               if(kmax1/=0)then
63                 kount=kount+1
64                 xp(kount)=x
65                 do i=1,nvar
66                   yp(i,kount)=y(i)
67                 end do
68               endif
69               return
70              endif
71             if(abs(hnext)<hmin) write(*,*) 'WARNING: stepsize smaller than minimum in odeint'
72             h=hnext
73            end do
74           write(*,*) 'WARNING: too many steps in odeint'
75           return
76           END
77