source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/oldies/twostepj.f @ 3947

Last change on this file since 3947 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 5.4 KB
Line 
1      SUBROUTINE INTEGRATE( TIN, TOUT )
2
3      INCLUDE 'KPP_ROOT_params.h'
4      INCLUDE 'KPP_ROOT_global.h'
5
6C TIN - Start Time
7      KPP_REAL TIN
8C TOUT - End Time
9      KPP_REAL TOUT
10      EXTERNAL ITER
11      KPP_REAL T
12      KPP_REAL V(NVAR), VOLD(NVAR), VNEW(NVAR)
13      KPP_REAL startdt, hmin, hmax, h 
14
15      INTEGER    INFO(5)
16
17      INFO(1) = Autonomous
18      h = hmin
19     
20c Number of Jacobi-Seidel iterations     
21      numit = 3   
22
23     
24      DO i=1,NVAR
25        RTOL(i) = 1.e-2
26      ENDDO     
27
28      CALL twostepj(NVAR,TIN,TOUT,h,hmin,hmax,
29     +                   VOLD,VAR,VNEW, 
30     +                   ATOL,RTOL,numit,
31     +                   nfcn,naccpt,nrejec,nstart,startdt,ITER)
32
33
34      RETURN
35      END
36     
37     
38
39      SUBROUTINE ITER(n,T,y,yp,yl)
40      INCLUDE 'KPP_ROOT_params.h'
41      INCLUDE 'KPP_ROOT_global.h'
42      REAL*8 T, y(NVAR), yp(NVAR), yl(NVAR)
43      TOLD = TIME
44      TIME = T
45      CALL Update_SUN()
46      CALL Update_RCONST()
47      CALL FunSPLIT_VAR(y,Rad,yp,yl)
48      TIME = TOLD
49      RETURN
50      END
51     
52
53      subroutine twostepj(n,t,te,dt,dtmin,dtmax,
54     +                   yold,y,ynew,
55     +                   atol,rtol,numit,
56     +                   nfcn,naccpt,nrejec,nstart,startdt,ITER)
57      implicit real*8 (a-h,o-z)
58      external ITER
59      integer    n,numit,nfcn,naccpt,nrejec,nstart,i,j
60      real*8     t,te,dt,dtmin,dtmax,startdt,ytol,
61     +           ratio,dtold,a1,a2,c,cp1,dtg,errlte,dy
62      real*8       yold(n),y(n),ynew(n),yp(n),yl(n),
63     +           work(n),sum(n),atol(n),rtol(n)
64      logical    accept,failer,restart
65
66c
67c     Initialization of counters, etc.
68
69      naccpt=0
70      nrejec=0
71      nfcn=0
72      nstart=0
73      failer=.false.
74      restart=.false.
75      accept=.true.
76
77c     Initial stepsize computation. 
78   
79   10 if (dtmin.eq.dtmax) then
80       nstart=1
81       dt=min(dtmin,(te-t)/2)
82       goto 28
83      endif
84      CALL ITER(n,t,y,yp,yl)
85      nfcn=nfcn+1
86      dt=te-t
87      do 20 i=1,n
88       ytol=atol(i)+rtol(i)*abs(y(i))
89       dy=yp(i)-y(i)*yl(i)
90       if (dy.ne.0.0) dt=min(dt,ytol/abs(dy))
91   20 continue
92   25 nstart=nstart+1
93      if (restart) dt=dt/10.0
94      restart=.true.
95      dt=max(dtmin,min(dt,dtmax))
96      CALL FIT(t,te,dt)
97      dt=min(dt,(te-t)/2)
98      startdt=dt
99
100c     The starting step is carried out, using the implicit Euler method.
101
102   28 do 30 i=1,n
103       ynew(i)=y(i)
104       yold(i)=y(i)
105       sum(i)=y(i)
106   30 continue
107      do 40 i=1,numit
108       CALL ITER(n,t+dt,ynew,yp,yl)
109       do i2i=1,n
110         ynew(i2i) = (sum(i2i) + dt*yp(i2i))/(1.+dt*yl(i2i))
111       end do
112       nfcn=nfcn+1
113   40 continue
114      naccpt=naccpt+1
115      t=t+dt
116      do 50 j=1,n
117       y(j)=ynew(j)
118   50 continue
119
120c     Subsequent steps are carried out with the two-step BDF method.
121
122      dtold=dt
123      ratio=1.0
124   60 continue
125      c=1.0/ratio
126      cp1=c+1.0
127      a1=((c+1.0)**2)/(c*c+2.0*c)
128      a2=-1.0/(c*c+2.0*c)
129      dtg=dt*(1.0+c)/(2.0+c)
130      do 70 j=1,n
131       sum(j)=a1*y(j)+a2*yold(j)
132       ynew(j)=max(0.0,y(j)+ratio*(y(j)-yold(j)))
133   70 continue
134      do 80 i=1,numit
135       CALL ITER(n,t+dt,ynew,yp,yl)
136       do i2i=1,n
137        ynew(i2i) = (sum(i2i) + dtg*yp(i2i))/(1.+dtg*yl(i2i))
138       end do
139       nfcn=nfcn+1
140   80 continue
141
142c     If stepsizes should remain equal, stepsize control is omitted.
143
144      if (dtmin.eq.dtmax) then
145       t=t+dtold
146       naccpt=naccpt+1
147       do 85 j=1,n
148        yold(j)=y(j)
149        y(j)=ynew(j)
150   85  continue
151       if (dt.ne.dtold) then
152        t=t-dtold+dt
153        goto 120
154       endif
155       dt=min(dtold,te-t)
156       ratio=dt/dtold
157       if (t.ge.te) goto 120
158       goto 60
159      endif
160
161c     Otherwise stepsize control is carried out.
162
163      errlte=0.0
164      do 90 i=1,n
165       ytol=atol(i)+rtol(i)*abs(y(i))
166       errlte=max(errlte,abs(c*ynew(i)-cp1*y(i)+yold(i))/ytol)
167   90 continue
168      errlte=2.0*errlte/(c+c*c)
169      CALL NEWDT(t,te,dt,dtold,ratio,errlte,accept,
170     +           dtmin,dtmax)
171
172c     Here the step has been accepted. 
173
174      if (accept) then
175 201   format(2(E24.16,1X))
176       failer=.false.
177       restart=.false.
178       t=t+dtold
179       naccpt=naccpt+1
180       do 100 j=1,n
181        yold(j)=y(j)
182        y(j)=ynew(j)
183  100  continue
184       if (t.ge.te) goto 120
185       goto 60
186      endif
187
188c     A restart check is carried out.
189     
190      if (failer) then
191       nrejec=nrejec+1
192       failer=.false.
193       naccpt=naccpt-1
194       t=t-dtold
195       do 110 j=1,n
196        y(j)=yold(j)
197  110  continue
198       goto 25
199      endif
200
201c     Here the step has been rejected.
202
203      nrejec=nrejec+1
204      failer=.true.
205      goto 60
206     
207c     End of TWOSTEP.
208  120 end
209c=====================================================================
210
211      subroutine NEWDT(t,te,dt,dtold,ratio,errlte,
212     +                 accept,dtmin,dtmax)
213      real*8    t,te,dt,dtold,ratio,errlte,ts,dtmin,dtmax
214      logical accept
215      if (errlte.gt.1.0.and.dt.gt.dtmin) then
216       accept=.false.
217       ts=t
218      else
219       accept=.true.
220       dtold=dt
221       ts=t+dtold
222      endif
223      dt=max(0.5,min(2.0,0.8/sqrt(errlte)))*dt
224      dt=max(dtmin,min(dt,dtmax))
225      CALL FIT(ts,te,dt)
226      ratio=dt/dtold
227      end
228
229      subroutine FIT(t,te,dt)
230      real*8    t,te,dt,rns
231      integer   ns
232      rns=(te-t)/dt
233      if (rns.gt.10.0) goto 10
234      ns=int(rns)+1
235      dt=(te-t)/ns
236      dt=(dt+t)-t
237   10 return
238      end
239
240
241
242C End of MAIN function                                             
243C **************************************************************** 
244
Note: See TracBrowser for help on using the repository browser.