Forgot to load lapack in a few examples
[maxima.git] / share / odepack / fortran / opkdemo6.f
blob2cfc20bdd6f7f0646008095fdbd3094afefdefbe
1 program opkdemo6
2 c-----------------------------------------------------------------------
3 c Demonstration program for the DLSODKR package.
4 c This is the version of 15 June 2001.
6 c This version is in double precision.
8 c An ODE system is generated from the following 2-species diurnal
9 c kinetics advection-diffusion PDE system in 2 space dimensions:
11 c dc(i)/dt = Kh*(d/dx)**2 c(i) + V*dc(i)/dx + (d/dz)(Kv(z)*dc(i)/dz)
12 c + Ri(c1,c2,t) for i = 1,2, where
13 c R1(c1,c2,t) = -q1*c1*c3 - q2*c1*c2 + 2*q3(t)*c3 + q4(t)*c2 ,
14 c R2(c1,c2,t) = q1*c1*c3 - q2*c1*c2 - q4(t)*c2 ,
15 c Kv(z) = Kv0*exp(z/5) ,
16 c Kh, V, Kv0, q1, q2, and c3 are constants, and q3(t) and q4(t)
17 c vary diurnally. The species are oxygen singlet and ozone.
18 c The problem is posed on the square
19 c 0 .le. x .le. 20, 30 .le. z .le. 50 (all in km),
20 c with homogeneous Neumann boundary conditions, and for time t in
21 c 0 .le. t .le. 86400 sec (1 day).
23 c The PDE system is treated by central differences on a uniform
24 c 10 x 10 mesh, with simple polynomial initial profiles.
26 c The problem is solved with DLSODKR, with the BDF/GMRES method and
27 c the block-diagonal part of the Jacobian as a left preconditioner.
28 c At intervals of 7200 sec (2 hrs), output includes sample values
29 c of c1, c2, and c2tot = total of all c2 values.
31 c Roots of the function g = d(c2tot)/dt are found, i.e. the points at
32 c which the total c2 (ozone) is stationary.
34 c Note: The preconditioner routines call LINPACK routines DGEFA, DGESL,
35 c and the BLAS routines DCOPY and DSCAL.
36 c-----------------------------------------------------------------------
37 external fdem, jacbd, solbd, gdem
38 integer mx, mz, mm, iout, istate, iwork, jacflg, jpre, jroot,
39 1 jx, jz, leniw, lenrw, liw, lrw, mf, neq, nst, nsfi, nfe,
40 2 nge, npe, njev, nps, nni, nli, ncfn, ncfl
41 double precision q1,q2,q3,q4, a3,a4, om, c3, dz, hdco, vdco, haco,
42 1 dkh, vel, dkv0, halfda, pi, twohr, rtol, floor,
43 2 dx, atol, t, tout, x, cx, z, cz, y, rwork, c2tot, avdim
44 common /pcom/ q1,q2,q3,q4,a3,a4,om,c3,dz,hdco,vdco,haco,mx,mz,mm
45 dimension y(2,10,10), rwork(4264), iwork(235)
46 dimension neq(1), rtol(1), atol(1), jroot(1)
47 data dkh/4.0d-6/, vel/0.001d0/, dkv0/1.0d-8/, halfda/4.32d4/,
48 1 pi/3.1415926535898d0/, twohr/7200.0d0/, rtol/1.0d-5/,
49 2 floor/100.0d0/, lrw/4264/, liw/235/, mf/22/, jpre/1/, jacflg/1/
51 c Load Common block of problem parameters.
52 mx = 10
53 mz = 10
54 mm = mx*mz
55 q1 = 1.63d-16
56 q2 = 4.66d-16
57 a3 = 22.62d0
58 a4 = 7.601d0
59 om = pi/halfda
60 c3 = 3.7d16
61 dx = 20.0d0/(mx - 1.0d0)
62 dz = 20.0d0/(mz - 1.0d0)
63 hdco = dkh/dx**2
64 haco = vel/(2.0d0*dx)
65 vdco = (1.0d0/dz**2)*dkv0
66 c Set other input arguments.
67 atol(1) = rtol(1)*floor
68 neq(1) = 2*mx*mz
69 iwork(1) = 8*mx*mz
70 iwork(2) = neq(1)
71 iwork(3) = jpre
72 iwork(4) = jacflg
73 t = 0.0d0
74 tout = 0.0d0
75 istate = 1
76 c Set initial profiles.
77 do 20 jz = 1,mz
78 z = 30.0d0 + (jz - 1.0d0)*dz
79 cz = (0.1d0*(z - 40.0d0))**2
80 cz = 1.0d0 - cz + 0.5d0*cz**2
81 do 10 jx = 1,mx
82 x = (jx - 1.0d0)*dx
83 cx = (0.1d0*(x - 10.0d0))**2
84 cx = 1.0d0 - cx + 0.5d0*cx**2
85 y(1,jx,jz) = 1.0d6*cx*cz
86 y(2,jx,jz) = 1.0d12*cx*cz
87 10 continue
88 20 continue
90 c Write heading, problem parameters, solution parameters.
91 write(6,30) mx, mz, mf, rtol(1), atol(1)
92 30 format('Demonstration program for DLSODKR package'//
93 1 '2D diurnal kinetics-transport PDE system with 2 species'/
94 2 'Spatial mesh is',i3,' by',i3/'Method flag is mf =',i3,
95 3 ' Tolerances are rtol =',d8.1,' atol =',d8.1/
96 4 'Left preconditioner uses block-diagonal part of Jacobian'/
97 5 'Root function finds stationary points of total ozone,'/
98 6 ' i.e. roots of (d/dt)(sum of c2 over all mesh points)'/)
100 c Loop over output points, call DLSODKR, print sample solution values.
101 do 70 iout = 1,13
102 40 call dlsodkr (fdem, neq, y, t, tout, 1, rtol, atol, 1, istate,
103 1 0, rwork, lrw, iwork, liw, jacbd, solbd, mf, gdem, 1, jroot)
104 write(6,50) t,iwork(11),iwork(14),rwork(11)
105 50 format(/' t =',d10.3,4x,'no. steps =',i5,
106 1 ' order =',i2,' stepsize =',d10.3)
107 call c2sum (y, mx, mz, c2tot)
108 write(6,60) y(1,1,1), y(1,5,5), y(1,10,10),
109 1 y(2,1,1), y(2,5,5), y(2,10,10)
110 60 format(' c1 (bot.left/middle/top rt.) =',3d12.3/
111 1 ' c2 (bot.left/middle/top rt.) =',3d12.3)
112 write(6,62)c2tot,jroot(1)
113 62 format(' total c2 =',d15.6,
114 1 ' jroot =',i2,' (1 = root found, 0 = no root)')
115 if (istate .lt. 0) then
116 write(6,65)istate
117 65 format('DLSODKR returned istate = ',i3)
118 go to 75
119 endif
120 if (istate .eq. 3) then
121 istate = 2
122 go to 40
123 endif
124 tout = tout + twohr
125 70 continue
127 c Print final statistics.
128 75 lenrw = iwork(17)
129 leniw = iwork(18)
130 nst = iwork(11)
131 nsfi = iwork(24)
132 nfe = iwork(12)
133 nge = iwork(10)
134 npe = iwork(13)
135 njev = iwork(25)
136 nps = iwork(21)
137 nni = iwork(19)
138 nli = iwork(20)
139 avdim = real(nli)/real(nni)
140 ncfn = iwork(22)
141 ncfl = iwork(23)
142 write (6,80) lenrw,leniw,nst,nsfi,nfe,nge,npe,njev,nps,nni,nli,
143 1 avdim,ncfn,ncfl
144 80 format(//' Final statistics:'/
145 1 ' rwork size =',i5,5x,' iwork size =',i4/
146 2 ' number of steps =',i5,5x,'no. fnal. iter. steps =',i5/
147 3 ' number of f evals. =',i5,5x,'number of g evals. =',i5/
148 4 ' number of prec. evals. =',i5,5x,'number of Jac. evals. =',i5/
149 5 ' number of prec. solves =',i5/
150 6 ' number of nonl. iters. =',i5,5x,'number of lin. iters. =',i5/
151 7 ' average Krylov subspace dimension (nli/nni) =',f8.4/
152 8 ' number of conv. failures: nonlinear =',i3,' linear =',i3)
153 c stop
156 subroutine fdem (neq, t, y, ydot)
157 integer neq, mx, mz, mm,
158 1 iblok, iblok0, idn, iup, ileft, iright, jx, jz
159 dimension neq(*)
160 double precision t, y(2,*), ydot(2,*),
161 1 q1, q2, q3, q4, a3, a4, om, c3, dz, hdco, vdco, haco,
162 2 c1, c1dn, c1up, c1lt, c1rt, c2, c2dn, c2up, c2lt, c2rt,
163 3 czdn, czup, horad1, hord1, horad2, hord2, qq1, qq2, qq3, qq4,
164 4 rkin1, rkin2, s, vertd1, vertd2, zdn, zup
165 common /pcom/ q1,q2,q3,q4,a3,a4,om,c3,dz,hdco,vdco,haco,mx,mz,mm
167 c Set diurnal rate coefficients.
168 s = sin(om*t)
169 if (s .gt. 0.0d0) then
170 q3 = exp(-a3/s)
171 q4 = exp(-a4/s)
172 else
173 q3 = 0.0d0
174 q4 = 0.0d0
175 endif
176 c Loop over all grid points.
177 do 20 jz = 1,mz
178 zdn = 30.0d0 + (jz - 1.5d0)*dz
179 zup = zdn + dz
180 czdn = vdco*exp(0.2d0*zdn)
181 czup = vdco*exp(0.2d0*zup)
182 iblok0 = (jz-1)*mx
183 idn = -mx
184 if (jz .eq. 1) idn = mx
185 iup = mx
186 if (jz .eq. mz) iup = -mx
187 do 10 jx = 1,mx
188 iblok = iblok0 + jx
189 c1 = y(1,iblok)
190 c2 = y(2,iblok)
191 c Set kinetic rate terms.
192 qq1 = q1*c1*c3
193 qq2 = q2*c1*c2
194 qq3 = q3*c3
195 qq4 = q4*c2
196 rkin1 = -qq1 - qq2 + 2.0d0*qq3 + qq4
197 rkin2 = qq1 - qq2 - qq4
198 c Set vertical diffusion terms.
199 c1dn = y(1,iblok+idn)
200 c2dn = y(2,iblok+idn)
201 c1up = y(1,iblok+iup)
202 c2up = y(2,iblok+iup)
203 vertd1 = czup*(c1up - c1) - czdn*(c1 - c1dn)
204 vertd2 = czup*(c2up - c2) - czdn*(c2 - c2dn)
205 c Set horizontal diffusion and advection terms.
206 ileft = -1
207 if (jx .eq. 1) ileft = 1
208 iright = 1
209 if (jx .eq. mx) iright = -1
210 c1lt = y(1,iblok+ileft)
211 c2lt = y(2,iblok+ileft)
212 c1rt = y(1,iblok+iright)
213 c2rt = y(2,iblok+iright)
214 hord1 = hdco*(c1rt - 2.0d0*c1 + c1lt)
215 hord2 = hdco*(c2rt - 2.0d0*c2 + c2lt)
216 horad1 = haco*(c1rt - c1lt)
217 horad2 = haco*(c2rt - c2lt)
218 c Load all terms into ydot.
219 ydot(1,iblok) = vertd1 + hord1 + horad1 + rkin1
220 ydot(2,iblok) = vertd2 + hord2 + horad2 + rkin2
221 10 continue
222 20 continue
223 return
226 subroutine gdem (neq, t, y, ng, gout)
227 integer neq, ng, mx, mz, mm,
228 1 iblok, iblok0, idn, iup, ileft, iright, jx, jz
229 dimension neq(*)
230 double precision t, y(2,*), gout(*),
231 1 q1, q2, q3, q4, a3, a4, om, c3, dz, hdco, vdco, haco,
232 2 c1, c2, c2dn, c2up, c2lt, c2rt, c2dot, czdn, czup, horad2,
233 3 hord2, qq1, qq2, qq4, rkin2, s, sum, vertd2, zdn, zup
234 common /pcom/ q1,q2,q3,q4,a3,a4,om,c3,dz,hdco,vdco,haco,mx,mz,mm
236 c This routine computes the rates for c2 and adds them.
238 c Set diurnal rate coefficient q4.
239 s = sin(om*t)
240 if (s .gt. 0.0d0) then
241 q4 = exp(-a4/s)
242 else
243 q4 = 0.0d0
244 endif
245 sum = 0.0d0
246 c Loop over all grid points.
247 do 20 jz = 1,mz
248 zdn = 30.0d0 + (jz - 1.5d0)*dz
249 zup = zdn + dz
250 czdn = vdco*exp(0.2d0*zdn)
251 czup = vdco*exp(0.2d0*zup)
252 iblok0 = (jz-1)*mx
253 idn = -mx
254 if (jz .eq. 1) idn = mx
255 iup = mx
256 if (jz .eq. mz) iup = -mx
257 do 10 jx = 1,mx
258 iblok = iblok0 + jx
259 c1 = y(1,iblok)
260 c2 = y(2,iblok)
261 c Set kinetic rate term for c2.
262 qq1 = q1*c1*c3
263 qq2 = q2*c1*c2
264 qq4 = q4*c2
265 rkin2 = qq1 - qq2 - qq4
266 c Set vertical diffusion terms for c2.
267 c2dn = y(2,iblok+idn)
268 c2up = y(2,iblok+iup)
269 vertd2 = czup*(c2up - c2) - czdn*(c2 - c2dn)
270 c Set horizontal diffusion and advection terms for c2.
271 ileft = -1
272 if (jx .eq. 1) ileft = 1
273 iright = 1
274 if (jx .eq. mx) iright = -1
275 c2lt = y(2,iblok+ileft)
276 c2rt = y(2,iblok+iright)
277 hord2 = hdco*(c2rt - 2.0d0*c2 + c2lt)
278 horad2 = haco*(c2rt - c2lt)
279 c Load all terms into c2dot and sum.
280 c2dot = vertd2 + hord2 + horad2 + rkin2
281 sum = sum + c2dot
282 10 continue
283 20 continue
284 gout(1) = sum
285 return
288 subroutine jacbd (f, neq, t, y, ysv, rewt, f0, f1, hl0, jok,
289 1 bd, ipbd, ier)
290 external f
291 integer neq, jok, ipbd(2,*), ier, mx, mz, mm,
292 1 iblok, iblok0, jx, jz, lenbd
293 dimension neq(*)
294 double precision t, y(2,*), ysv(*), rewt(*), f0(*), f1(*),
295 1 hl0, bd(2,2,*),
296 2 q1, q2, q3, q4, a3, a4, om, c3, dz, hdco, vdco, haco,
297 3 c1, c2, czdn, czup, diag, temp, zdn, zup
298 common /pcom/ q1,q2,q3,q4,a3,a4,om,c3,dz,hdco,vdco,haco,mx,mz,mm
300 lenbd = 4*mm
301 c If jok = 1, copy saved block-diagonal approximate Jacobian into bd.
302 if (jok .eq. 1) then
303 call dcopy (lenbd, bd(1,1,mm+1), 1, bd, 1)
304 go to 30
305 endif
307 c If jok = -1, compute and save diagonal Jacobian blocks
308 c (using q3 and q4 values computed on last f call).
309 do 20 jz = 1,mz
310 zdn = 30.0d0 + (jz - 1.5d0)*dz
311 zup = zdn + dz
312 czdn = vdco*exp(0.2d0*zdn)
313 czup = vdco*exp(0.2d0*zup)
314 diag = -(czdn + czup + 2.0d0*hdco)
315 iblok0 = (jz-1)*mx
316 do 10 jx = 1,mx
317 iblok = iblok0 + jx
318 c1 = y(1,iblok)
319 c2 = y(2,iblok)
320 bd(1,1,iblok) = (-q1*c3 - q2*c2) + diag
321 bd(1,2,iblok) = -q2*c1 + q4
322 bd(2,1,iblok) = q1*c3 - q2*c2
323 bd(2,2,iblok) = (-q2*c1 - q4) + diag
324 10 continue
325 20 continue
326 call dcopy (lenbd, bd, 1, bd(1,1,mm+1), 1)
327 c Scale by -hl0, add identity matrix and LU-decompose blocks.
328 30 temp = -hl0
329 call dscal (lenbd, temp, bd, 1)
330 do 40 iblok = 1,mm
331 bd(1,1,iblok) = bd(1,1,iblok) + 1.0d0
332 bd(2,2,iblok) = bd(2,2,iblok) + 1.0d0
333 call dgefa (bd(1,1,iblok), 2, 2, ipbd(1,iblok), ier)
334 if (ier .ne. 0) return
335 40 continue
336 return
339 subroutine solbd (neq, t, y, f0, wk, hl0, bd, ipbd, v, lr, ier)
340 integer neq, ipbd(2,*), lr, ier, mx, mz, mm, i
341 double precision t, y(*), f0(*), wk(*), hl0, bd(2,2,*),
342 2 v(2,*), q1, q2, q3, q4, a3, a4, om, c3, dz, hdco, vdco, haco
343 dimension neq(*)
344 common /pcom/ q1,q2,q3,q4,a3,a4,om,c3,dz,hdco,vdco,haco,mx,mz,mm
345 c Solve the block-diagonal system Px = v using LU factors stored in bd
346 c and pivot data in ipbd, and return the solution in v.
347 ier = 0
348 do 10 i = 1,mm
349 call dgesl (bd(1,1,i), 2, 2, ipbd(1,i), v(1,i), 0)
350 10 continue
351 return
354 subroutine c2sum (y, mx, mz, c2tot)
355 integer mx, mz, jx, jz
356 double precision y(2,mx,mz), c2tot, sum
357 c Sum the c2 values.
358 sum = 0.0d0
359 do 20 jz = 1,mz
360 do 20 jx = 1,mx
361 20 sum = sum + y(2,jx,jz)
362 c2tot = sum
363 return