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
19
20 =0.0
21 dxsav=0.0
22
23
24
25
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