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.
61 dx
= 20.0d0
/(mx
- 1.0d0
)
62 dz
= 20.0d0
/(mz
- 1.0d0
)
65 vdco
= (1.0d0
/dz**2
)*dkv0
66 c Set other input arguments.
67 atol
(1) = rtol
(1)*floor
76 c Set initial profiles.
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
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
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.
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
117 65 format('DLSODKR returned istate = ',i3
)
120 if (istate
.eq
. 3) then
127 c Print final statistics.
139 avdim
= real(nli
)/real(nni
)
142 write (6,80) lenrw
,leniw
,nst
,nsfi
,nfe
,nge
,npe
,njev
,nps
,nni
,nli
,
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
)
156 subroutine fdem
(neq
, t
, y
, ydot
)
157 integer neq
, mx
, mz
, mm
,
158 1 iblok
, iblok0
, idn
, iup
, ileft
, iright
, jx
, jz
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.
169 if (s
.gt
. 0.0d0
) then
176 c Loop over all grid points.
178 zdn
= 30.0d0
+ (jz
- 1.5d0
)*dz
180 czdn
= vdco*exp
(0.2d0*zdn
)
181 czup
= vdco*exp
(0.2d0*zup
)
184 if (jz
.eq
. 1) idn
= mx
186 if (jz
.eq
. mz
) iup
= -mx
191 c Set kinetic rate terms.
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.
207 if (jx
.eq
. 1) ileft
= 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
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
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.
240 if (s
.gt
. 0.0d0
) then
246 c Loop over all grid points.
248 zdn
= 30.0d0
+ (jz
- 1.5d0
)*dz
250 czdn
= vdco*exp
(0.2d0*zdn
)
251 czup
= vdco*exp
(0.2d0*zup
)
254 if (jz
.eq
. 1) idn
= mx
256 if (jz
.eq
. mz
) iup
= -mx
261 c Set kinetic rate term for 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.
272 if (jx
.eq
. 1) ileft
= 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
288 subroutine jacbd
(f
, neq
, t
, y
, ysv
, rewt
, f0
, f1
, hl0
, jok
,
291 integer neq
, jok
, ipbd
(2,*), ier
, mx
, mz
, mm
,
292 1 iblok
, iblok0
, jx
, jz
, lenbd
294 double precision t
, y
(2,*), ysv
(*), rewt
(*), f0
(*), f1
(*),
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
301 c If jok = 1, copy saved block-diagonal approximate Jacobian into bd.
303 call dcopy
(lenbd
, bd
(1,1,mm
+1), 1, bd
, 1)
307 c If jok = -1, compute and save diagonal Jacobian blocks
308 c (using q3 and q4 values computed on last f call).
310 zdn
= 30.0d0
+ (jz
- 1.5d0
)*dz
312 czdn
= vdco*exp
(0.2d0*zdn
)
313 czup
= vdco*exp
(0.2d0*zup
)
314 diag
= -(czdn
+ czup
+ 2.0d0*hdco
)
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
326 call dcopy
(lenbd
, bd
, 1, bd
(1,1,mm
+1), 1)
327 c Scale by -hl0, add identity matrix and LU-decompose blocks.
329 call dscal
(lenbd
, temp
, bd
, 1)
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
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
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.
349 call dgesl
(bd
(1,1,i
), 2, 2, ipbd
(1,i
), v
(1,i
), 0)
354 subroutine c2sum
(y
, mx
, mz
, c2tot
)
355 integer mx
, mz
, jx
, jz
356 double precision y
(2,mx
,mz
), c2tot
, sum
361 20 sum
= sum
+ y
(2,jx
,jz
)